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1 .  INTRODUCTION 


In  support  of  Space  Division's  concern  about  the  safety  of  orbiting  pay- 
loads,  the  Space  Hazards  Section  of  the  Astrodynamics  Department  was  tasked  to 
develop  an  analytical  model  for  debris  analysis  after  an  orbital  breakup.  The 
results  of  the  model  were  to  be  used  in  developing  a  computer  program  which 
could  examine  the  collision  hazard  to  any  spacecraft  from  the  cloud  of  parti¬ 
cles  resulting  from  an  orbital  breakup.  The  study  activity  consisted  of 
defining  the  requirements,  examining  the  available  hypervelocity  impact  data, 
and  developing  the  appropriate  breakup  models  which  could  be  used  to  determine 
the  fragment  population  and  velocity  distributions  for  particles  in  orbiting 
debris  clouds. 

Section  2  and  Appendix  A,  written  by  V.A.  Chobotov,  derive  the  analytical 
model  for  the  debris  cloud,  assuming  a  breakup  or  a  collision  with  a  space 
object  in  a  circular  orbit.  Linearized  equations  for  relative  motion  are  used 
to  determine  the  shape  and  volume  of  the  debris  cloud  for  an  initially  iso¬ 
tropic  distribution  of  particle  spread  velocities.  Spatial  density  is  obtained 
for  two  representative  breakup  models  and  the  collision  probability  of  a  parti¬ 
cle  and  a  resident  space  object  determined.  The  effects  of  earth's  oblateness 
on  the  long-term  evolution  of  the  cloud  are  examined. 

Section  3,  written  by  D.B.  Spencer,  compares  the  exact  with  the  linear 
approximation  results  from  Section  2.  It  is  shown  that  the  linear  approxima¬ 
tion  solution  generally  is  valid  for  the  case  of  low  particle  spread  velocities 
(<  100  ft/sec).  For  greater  velocities  (>  1000  ft/sec),  the  radial  and 
tangential  position  components  for  orbital  plane  ejections  begin  to  deteriorate 
as  time  increases.  The  differences  between  the  exact  and  linear  solutions  for 
particle  trajectories  are  illustrated. 

In  Section  A,  also  written  by  D.B.  Spencer,  the  volume  of  the  debris 
cloud  is  reexamined  with  the  inclusion  of  earth's  oblateness  and  atmospheric 
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drag  effects.  Time-dependent  functions  are  derived  which  model  the  changes  in 
the  cloud  profile  as  the  result  of  such  perturbations.  These  functions  are 
inputs  also  to  the  "DEBRIS"  computer  simulation  which  was  developed  at  The 
Aerospace  Corporation  for  the  collision  hazard  assessment  purpose. 

Section  5  and  Appendix  B,  contributed  by  D.L.  Schmitt,  describe  a  method 
for  determining  the  masses  and  velocities  of  fragments  resulting  from  a  hyper¬ 
velocity  collision  in  orbit  and  compute  their  orbital  parameters.  Tabular 
data  are  given  showing  the  percentages  of  fragments  that  reenter  or  stay  in 
orbit.  Plots  of  orbital  distributions  and  reentry  footprints  are  illustrated. 

Section  6,  written  by  R.P.  Gupta,  generalizes  the  fragment  mass  and 
velocity  computation  described  in  Section  5.  The  methodology  presented  in 
Section  6  can  be  used  to  determine  number,  mass,  size,  and  velocity  distribu¬ 
tion  of  fragments  for  different  mass  ratios  of  impacting  objects  and  includes 
the  effects  of  energy  loss  due  to  heat  and  light  generated  by  impact. 

Section  7,  authored  by  D.T.  Knapp,  develops  a  new  and  more  flexible  model 
for  spacecraft  collisions.  A  system  of  transfer  functions  characterizes  the 
collision  in  terms  of  the  incident  object  masses  and  velocities  and  a  set  of 
parameters  defined  to  reflect  the  several  degrees  of  freedom  of  the  system 
constrained  by  conservation  of  mass,  momentum,  and  energy.  This  allows 
sampling  of  many  parametric  values  for  a  statistical  or  sensitivity  analysis 
of  the  system.  Referred  to  as  the  Kinematic  Model,  it  overcomes  significant 
limitations  of  earlier  models.  It  has  been  used  by  The  Aerospace  Corporation 
to  model  tests  involving  on-orbit  collisions,  and  its  results  have  been 
verified  by  test  data. 

In  Section  8,  contributed  by  R.G.  Hopkins,  the  description  of  the  program 
DEBRIS  is  provided.  The  program  determines  the  intervals  during  which  a  space¬ 
craft  travels  through  an  expanding  cloud  and  calculates  the  probability  of 
collision  associated  with  each  transit. 
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Finally,  Section  9,  written  by  D.B.  Spencer,  examines  the  effects  of 
orbital  eccentricity  on  the  volume  of  the  debris  cloud.  Small  values  of 
eccentricity  are  added  to  the  originally  circular  orbit  of  the  disintegrating 
body  (as  was  assumed  in  Section  2);  by  using  a  differential  correction 
process,  the  changes  in  the  cloud  volume  are  determined. 

In  summary,  the  results  presented  in  this  report  represent  the  theoreti 
cal  background  and  description  of  the  DEBRIS  program  development  which  can  be 
used  to  assess  the  collision  hazard  for  resident  space  objects  following  a 
breakup  or  a  collision  of  an  object  in  orbit. 


2.  DYNAMICS  OF  DEBRIS  AND  THE  CONSEQUENCES  OF  ON-ORBIT  BREAKUPS 


2.1  INTRODUCTION 

Continuous  use  of  space  over  the  last  30  years  has  built  up  a  large  num¬ 
ber  of  objects  in  orbit,  the  majority  of  which  were  generated  by  explosions  of 
spacecraft  or  rocket  stages.  Late  1970s  studies  at  the  Johnson  Space  Center 
concluded  that  fragments  from  collisions  between  space  objects  would  be  a  major 
source  of  debris  (Ref.  1).  Studies  at  The  Aerospace  Corporation  examined  the 
collision  hazard  to  operational  spacecraft  from  space  debris  including  the 
effects  of  position  uncertainty  on  the  probability  of  collision  between  any  two 
objects  in  orbit  (Ref.  2).  Other  studies,  described  in  Reference  3,  considered 
the  distribution  of  some  5000  NORAD  Catalog  objects  as  a  function  of  altitude 
and  orbital  inclination.  Encounter  parameters  such  as  miss  distance  and  rela¬ 
tive  velocity  were  examined  by  computer  simulation  for  low-altitude  and  geosyn¬ 
chronous  orbit  spacecraft.  Representative  space  shuttle  and  geosynchronous 
mission  collision  probabilities  were  determined. 

One  of  the  early  studies  which  addressed  the  evolution  of  a  fragment 
cloud  in  orbit  is  discussed  in  Reference  4.  In  that  study,  the  time  and  place 
of  a  satellite  disintegration  were  determined  from  the  orbits  of  the  individual 
particles  (fragments)  obtained  by  observation.  Methods  of  statistical  mecha¬ 
nics  also  were  used  to  study  the  evolution  of  the  fragment  cloud  by  treating 
the  fragments  as  an  ensemble  of  noninteracting  particles.  The  spatial  density 
was  calculated  as  a  function  of  position,  time,  and  initial  velocity 
distribution. 

This  study  considers  the  problem  of  debris  cloud  evolution  by  examining 
representative  particle  trajectories.  Linearized  equations  for  relative  motion 
in  orbit  are  used  to  obtain  the  trajectories  of  particles  with  specified 
initial  velocity  distributions  in  three  orthogonal  planes.  The  volume  of  the 
cloud  is  computed  analytically  as  a  function  of  time,  and  the  spatial  density 
is  calculated  for  representative  breakup  models.  Long-term  effects  due  to 
earth's  oblateness  are  evaluated,  and  the  near-  and  long-term  collision 
hazards  for  representative  spacecraft  are  examined. 
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2.2 


ANALYSIS 


Consider  an  explosion  or  a  collision  event  in  a  circular  orbit  such  as 
that  illustrated  in  Figure  1.  An  orbiting  orthogonal  reference  fraune  xyz  is 
centered  at  the  origin  of  the  event  at  time  t  =  0  such  that  x  is  directed 
opposite  to  the  orbital  velocity  vector,  y  is  directed  along  the  outward 
radius,  and  z  completes  the  triad  (along  the  normal  to  the  orbit  plane).  The 
linearized  rendezvous  equations  (Ref.  5)  can  be  used  to  determine  the  position 
of  a  particle  leaving  the  origin  of  the  coordinate  frame  with  a  velocity  Av; 
they  are  of  the  form 


X  =  +  —  sin  0)x  +  -(1  -  cos  0)y 

\  u  u  /  o  M  c 


y  =  -(cos  0  -  l)x  +  -^  sin  0 

■'  bi  Ob) 


z 

_c 

(l) 


o  •  a 
z  =  —  Sin  0 


(1) 


V.  A  ^*2  *2  *2.1/2 

where  Av  =  (x  +  y  +  z  ) 
o  "^o  0 


Figure  1.  Cloud  Dynamics 
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The  x,y,z  coordinates  represent  particle  position  at  time  t  =  9/w, 

where  0  is  the  in-orbit  plane  angle  and  a  is  the  angular  rate  of  the  circular 

orbit.  The  x  ,  y  ,  z  terms  are  initial  velocity  components  imparted  to  the 
o  o  o 

particle  along  the  x,y,z  axes,  respectively.  It  is  assumed  that  Av  «  v. 


Equation  (1)  can  be  normalized  with  respect  to  Av/u  as  follows: 


^  =  (-30  +  U  sin  0)h  +  2(1  -  cos  0)r 
Av 


=  X 


^  =  2(cos  0  -  l)h  +  (sin  0)r 
Av 

=  Y 


Z(i» 

Av 


n  sin  0 


»  Z 


> 


(2) 


where 


•  •  •  2  2  2 

h  =  X  /Av,  ray  /Av,  n  =  z  /Av,  and  h  r  +  n  *  1  (3) 

o  0  o 

In  matrix  form 


X 

ail  *12  *13 

h 

Y 

s 

*21  *22  *23 

r 

z 

_  *31  *32  *33. 

n 

=  [M] 


h 

r 

n 


(4) 


where 


a 


11  " 
®21  “ 
®31  = 


(-3  0  +  A  sin  0) 
2(cos  0  -1) 

0 


a^2  “  2(1  -  cos  0) 

a22  “  sin  ® 

*32  •  “ 


*13  ■  ° 

*23  '  “ 
a^j  =  sin  0 


Eqviations  (4)  can  be  plotted  as  a  function  of  0  for  different  values 
of  h,  r,  and  n  satisfying  condition  (3).  If,  for  example,  the  initial  velocity 
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Av  distribution  for  the  particles  is  circular  in  the  x,y  plane,  then  h  =  cos 

2  2 

r  =  sin  A<^,  0  <  A(^  <  360®,  and  h  +  r  =1  with  n  =  0.  The  resultant  cloud 
outline  is  illustrated  in  Figure  2  for  several  values  of  0. 


y 


Figure  2.  Debris  Cloud  in  Orbit  Plane 

Representative  cloud  outlines  in  nondimensional  units  are  shown  in 
Figure  3.  This  figure  illustrates  the  cloud  outlines  in  the  plane  of  the 
orbit  (xy)  at  four  different  times  after  breakup.  The  straight  line  configura¬ 
tion  results  after  one  revolution  (period  of  the  nominal  circular  orbit).  It 
shows  that  an  equal  number  of  particles  are  leading  and  lagging  as  expected 
from  the  uniform  velocity  distribution  assumed  initially  at  6  =  0®. 
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Figure  3.  Cloud  Contours  in  Orbit  Plane  (to  scale) 

For  a  velocity  distribution  that  is  circular  in  the  xz  plane,  the  cloud 
outlines  in  Figure  4  represent  the  contours  in  the  cross-track  (xz)  plane.  A 
straight  line  configuration  occurs  at  0  »  180*  and  360*  where  all  particles 
cross  the  orbit  plane. 
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Figure  4.  Cloud  Contours  in  Cross-Track  Plane  (to  scale) 
2 . 3  VOLUME  OF  DEBRIS  CLOUD 


The  volume  of  the  debris  cloud  can  be  expressed  analytically  in  terms 
of  the  elements  of  the  transition  matrix  M  in  Eq.  (4).  If,  for  example,  h,r,n 
are  the  orthogonal  components  of  a  sphere  of  unit  radius,  then  the  determinant 
of  M  represents  the  volume  of  the  cloud  at  any  time  t  =  0/u.  The  positive 
cloud  volume,  normalized  to  (Av/m)  ,  is  of  the  form 


VOLUME  =  ~  |det  [M] | 

=  |(ad  +  b^)dl 


(5) 
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where  a  =  b  =  aj^2»  ^  -  ®22’  where  Uk/3  is  the  volume  of  a  unit 

sphere.  A  further  discussion  of  this  theory  is  presented  in  Appendix  A. 
Equation  (5)  is  plotted  in  Figure  5  where  a  linear  approximation  for  the 
volume  also  is  shown. 

In  units  of  (Av/u),  the  volume  becomes 

.  3 

VOL  =  VOLUME  (— )  (6) 

u 


Figure  5.  Cloud  Volume  versus  Time 

For  a  circular,  low-altitude  orbit  with  Av  =  100  and  200  m/s  and 
-3  -1 

(I)  =  1.1  X  10  8  ,  the  volume  of  the  debris  cloud  is  illustrated  in  Figure  6. 

Note  that  the  volume  vanishes  at  the  integral  values  of  9  s  2^  and  6  »  x  and 
for  9  between  500®  and  550®  and  again  near  9  =  900®.  The  latter  zeros  are 
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caused  by  the  area  of  the  cloud  in  the  orbit  (xy)  plane  collapsing  to  a  line 
due  to  the  linearization  of  the  equations  of  motion.  A  condition  that  the 
area  not  vanish  anywhere  except  at  0  =  0,  ir,  2ir,  et  cetera,  can  be  satisfied  by 

ad  +  b^  >  0  (7) 

or 

lad|  +  b^  >  0  (8) 

which  results  in  the  cloud  volume  function  shown  in  Figure  7  generated  using 
Eq.  (5)  with  ad  >  0,  which  ensures  that  the  condition  (8)  is  always  satis¬ 
fied.  This  approximation  for  the  cloud  volume  improves  as  0  increases  when 

|ad|  »  b^  (9) 


2. A  SPATIAL  DENSITY 


Assuming  uniform  distribution  of  the  particles  in  the  cloud,  the  number 
of  which  is  to  be  determined  later,  the  density  is  of  the  form 

p  =  N/VOL  (10) 

where  N  is  the  number  of  particles  in  the  cloud.  The  accuracy  of  this  result 
is  reasonably  good  for  low  values  of  the  particle  spread  velocities  (e.g.,  Av 
<  100  m/s). 

A  "mean"  value  for  p  may,  for  example,  be  obtained  approximately  as 

^av  ~  (VOL) 

av 

where 

(VOL)  . 

av  6.65  u 
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THETA  (deg) 


Figure  7.  Cloud  Volume  versus  Time  (corrected) 


r 

I 


is  an  arbitrarily  assumed  linear  function  of  0,  as  can  be  seen  from  Figure  5. 
The  volume  of  the  initial  spherical  cloud  in  this  case  increases  linearly  with 
time.  Such  an  approximation  is  valid  up  to  a  quarter  revolution  (9  =  90®) 
when  the  "mean"  volume  is 


(VOL)  =  13.53(— )  (12) 

av  u 

-3  -1 

Thus,  for  example,  if  £iv  =  100  m/s,  to  =  1.1  x  10  s  ,  then 

(VOL)  =  1.016  X  lO^(km)^  (13) 

av 

3 

For  Av  =  200  m/s,  the  volume  is  a  factor  of  (2)  greater.  For 
Av  =  100  m/s,  the  spherical  cloud  diameter  at  1/4  revolution  (0  =  90®)  is  269  km 
which  compares  with  285  km  if  obtained  as  a  linear  function  of  Av. 

2.5  SHORT-TERM  COLLISION  HAZARD 


The  probability  that  a  spacecraft  will  collide  with  a  fragment  while 

passing  through  a  debris  cloud  is  proportional  to  the  spatial  density  in  the 

cloud,  spacecraft  projected  area.  A,  spacecraft  velocity  relative  to 

the  cloud,  V„,  and  the  time,  t,  spent  in  the  cloud.  Thus 
K 

p(col)  =  p^^AVj^t  (14) 

where  p(col)  =  collision  probability  per  pass. 

The  product  Vj^t  is  the  path  length  through  the  cloud.  It  is  the 
diameter  D  of  a  spherical  cloud  for  a  spacecraft  passing  through  the  center  of 
the  cloud.  The  probability  of  collision,  then,  is  of  the  form 

p(col)/A  =  p  D  (15) 

fl  V 
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where 


6 (VOL) 

_ av 

ir 

Evaluation  of  Eq.  (15)  requires  knowledge  of  N  and  Av  to  obtain 
and  (VOL)^^.  Laboratory  hypervelocity  impact  experiments,  such  as 
are  described  in  Reference  6,  for  example,  have  shown  that  the  number  versus 
size  distribution  for  ejecta  fragments  is  of  the  form 

N  .  0.4478(^ 
e 

where  N  is  the  cumulative  number  of  ejecta  with  mass  M  or  greater  and 
2 

M  =  M  v  where  M  is  the  mass  of  the  smaller  body  (projectile)  and  v  is  the 
e  p  p 

collision  velocity. 

Equation  (16)  and  Figure  8  from  Reference  6,  which  illustrate  fragment 
velocity  distribution  from  one  laboratory  test  with  an  impact  velocity  of 
3.5  kffl/s,  were  used  to  obtain  the  distributions  of  fragments  in  Table  1  from 
an  assumed  collision  of  two  objects  in  orbit.  Each  of  the  three  clouds  corres¬ 
ponds  to  a  different  particle  size. 

A  second  distribution  is  illustrated  in  Table  2  which  specifies  a  range 
of  particle  sizes  for  each  spread  velocity  group. 

The  probability  of  collision  of  a  spacecraft  and  a  debris  particle  based 
on  Eq.  (15)  is  shown  in  Figure  9.  The  curves  labeled  NASA  and  Aerospace 
correspond  to  the  distributions  of  debris  particles  in  Tables  1  and  2,  respec¬ 
tively. 


(16) 


The  results  show  that  the  collision  hazard  decreases  rapidly  after  the 
event  (0  =  0  in  Fig.  9),  but  that  the  magnitude  of  the  hazard  is  greatly 
dependent  on  the  assumed  distribution.  The  probabilities  for  each  of  the 
three  cloud  distributions  in  Tables  1  and  2  were  added  to  obtain  the  results 
in  Figure  9. 
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Figure  8.  Velocity  Distribution  (Ref.  6) 


Table  1.  NASA — ^Velocity  Distribution 


Cloud 

Number  of  Particles 
(N) 

Spread  Velocity 
Av(m/s) 

Particle  Size 
d(cm) 

1 

200 

20 

10 

2 

20000 

200 

1 

3 

3000000 

1000 

0.1 

Table  2.  The  Aerospace  Corporation — Velocity  Distribution  (Kinematic  model) 


Cloud 

Number  of  Particles 
(N) 

Spread  Velocity 
Av(m/s) 

Particle  Size 
d(cm) 

1 

200000 

0  to  100 

0.1  to  120 

2 

500000 

0  to  800 

0.1  to  60 

3 

1000000 

0  to  2000 

0.1  to  6 
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Figure  9.  Probability  of  Collision  per  Pass 
2.6  LONG-TERM  EFFECTS 

The  debris  particles  in  the  cloud  tend  to  spread  along  the  circumfer¬ 
ence  of  the  nominal  orbit,  in  time  assuming  the  shape  of  a  torus  with  two 
"pinch"  points  as  illustrated  in  Figure  10.  The  pinch  points  result  from  the 
orbital  intersections  of  the  debris  particles  where,  theoretically,  the  volume 
of  the  cloud  becomes  zero.  All  particles  pass  through  the  point  of  disinte¬ 
gration  (0  =  0)  at  different  times.  They  also  pass  through  the  orbit  plane 
at  9  =  180®  along  a  line  on  the  radius  vector.  The  probability  of  collision 
near  the  pinch  points  can  be  much  greater  where  the  cloud  volume  initially  is 
very  small.  The  effect  of  orbit  perturbations  (e.g.,  earth's  oblateness,  air 
drag,  etc.)  is  to  slowly  increase  the  volume  at  the  pinch  points  and  thus  de¬ 
crease  the  probability  of  collision.  In  the  long  term,  the  motion  of  the  line 
of  apsides  and  nodal  drift  tend  to  widen  the  pinch  points  and  spread  the  cloud 
envelope  until  it  completely  envelops  the  earth.  In  this  steady-state  condi¬ 
tion,  the  collision  hazard  is  reduced  to  a  minimum  and  can  be  compared  with 
the  existing  background  environment,  i.e.,  micrometeoroids,  man-made  debris, 
and  so  forth.  An  evaluation  of  the  apsidal  and  nodal  drifts  is  considered  in 
the  following  sections. 


If  the  particles'  individual  orbital  parameters  are  unknown,  the  expected 
mean  values  for  these  parameters  can  be  obtained  probabilistically  for  any 
specified  spread  velocity  distribution.  A  change  in  the  particle  semi-major 
axis  Aa  can  be  computed  as  a  function  of  the  spread  velocity  vector,  Av,  and 
the  orbital  velocity  v  geometry.  This  can  be  done  from  a  functional  relation¬ 
ship  between  Aa,  Av,  and  the  angle  between  v  and  Av,  as  follows. 


Consider  the  specific  energy  of  a  particle  orbit  of  the  form 


E  = 


y_ 

2a 


(17) 


where  yi  =  gravitational  constant.  Taking  differentials 


2a 


(18) 
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and 


AE  _  -Aa 
E  “  a 

This  equation  requires  an  explicit  expression  for  AE  which  is  the  change  in 
the  particle  orbit  energy  as  a  result  of  the  change  in  its  velocity  Av. 

This  is  of  the  form 

AE  =  |(v  +  Av)^  -  I  (v)^ 

=  |(v  +  Av)  •  (v  +  Av)  -  — 


=  V  Av  cos  ^  + 


»  V  Av  cos  4  for  —  «  1 

^  V 


Thus,  from  Eqs.  (17)  and  (18) 


V  Av  cos 
E 


for  ^  «  1 
V 


2av  Av  cos 
-VI 


Aa  av  Av  cos  4 
2a  ~  p 

„  Av 

»  —  cos  4 

V 

,  /  ^l/2 

since  V  =  (p/a) 

A  general  approximate  relationship  between  Aa,  Av,  and  4 
small  values  of  8v,  therefore,  is  of  the  form 


Sa  s  $v  cos  4 
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where 


&a 


Aa 

2a 


Equation  (23)  describes  the  internal  structure  of  the  cloud.  For  an 

2 

exact  relationship,  the  (Av)  /2  term  in  Eq.  (20)  must  be  retained. 


2.8  PROBABILITY  DISTRIBUTION 


The  functional  relationship  between  Av,  v,  and  4*  expressed  in  Eq.  (23) 
is  plotted  in  Figure  11.  Equation  (23)  or  Figure  11  can  be  used  to  obtain 


V 


Figure  11.  Internal  Structure  of  Cloud 
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orbital  parameter  probability  distributions  for  groups  of  particles  with 
specified  spread  velocity  ranges.  If,  for  example,  an  isotropic  spread 
velocity  distribution  is  assumed  as  shown  in  Figure  12,  then  the  probability 
that  Av  will  be  within  an  angle  of  v  is 


2A„ 


P  = 


SPH 


4-tr  V  (1  -  cos 
2 

4ir  V 

=  1  -  cos  ^ 


(24) 


where  A„  is  the  area  of  a  spherical  zone  defined  by  the  cone  angle  and  A 

I*  brn 

is  the  area  of  the  sphere  of  radius  v. 


y 


Figure  12.  Isotropic  Spread  Velocity  Distribution 
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Using  Eqs.  (23)  and  (2A),  we  find  that  the  probability  distribution 
function  is  of  the  form 


p  ~  I 


Sa 

Sv 


(25) 


or 


6a  =  (1  -  p)6v 


(26) 


Thus,  the  mean  or  the  "expected"  change  (corresponding  to  p  =  0.5)  in 

Sa  is 

Sa  =  i  (27) 

where  +  corresponds  co  whether  ^  is  less  or  greater  than  tr/Z.  One-half  of 
the  particle  orbits,  therefore,  will  have  a  change  in  the  semi-major  axes  of 


<Aa>  a  +  a(— ) 

-  V 


(28) 


Equation  (25)  is  plotted  in  Figure  13  with  Sv  as  a  parameter. 


Sa 


Figure  13.  Probability  Distribution  for  Sa  =  Aa/2a 
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EARTH'S  OBLATENESS  EFFECTS 


The  primary,  long-term  orbit  perturbation  effect  on  the  evolution  of 
the  debris  cloud  is  earth's  oblateness.  The  term  in  the  earth's  gravita¬ 
tional  potential  accounts  for  earth's  oblateness  and  causes  apsidal  and  nodal 
rotation.  This  in  turn  results  in  a  constantly  increasing  volume  of  the  cloud. 
In  order  to  determine  the  cloud  growth  rate  due  to  it  is  necessary  to 
determine  the  inclination  i,  semi-major  axis  a,  and  eccentricity  e  for  each 
particle  orbit.  The  apsidal  rotation  w  and  nodal  regression  5  are  then 
given  by  the  equations 


where 


(d  =  K(2  -  ^  sin^  i) 

(29) 

^  =  -K  cos  i 

(30) 

9.9639  (\\  . 

earth  eqtiatorial  radius 

(31) 

The  initial,  pinched  ring  shape  of  the  cloud  (such  as  is  shown  in 
Fig.  10)  eventually  spreads  arotmd  the  earth  due  to  particle  orbit  apsidal  and 
nodal  rotation  effects.  The  rate  at  which  this  steady-state  condition  is 
approached  can  be  obtained  from  the  mean  expected  values  <«>,  <§>  for  the 
apsidal  and  nodal  rotation  rates,  respectively. 


Consider,  for  example,  the  case  of  isotropic  spread  velocity  distribu¬ 
tion  in  a  low-altitude  circular  orbit  with  the  following  parameters: 


Av  s  100  m/s 
a  =  6924  km 

V  =  7.63  km/s  (32) 

i  =  98* 
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For  this  case,  5v  =  Av/v  =  0.0131,  and  the  probability  density 
distribution  for  5a  is  as  shown  by  the  dashed  line  in  Figure  14  (see  Section 
3).  From  Figure  14  or  Eq.  (28) 

<da>  =  +  a  5v 

=  +  90.7  km  (33) 

Particles  with  the  positive  mean  change  of  the  semi-major  axis  (one- 
half  of  all  fragments)  have  orbits  with  higher  energy  than  that  of  the  parent 
body,  while  those  with  negative  have  lower  energy.  The  corresponding  mean 
apsidal  and  nodal  rotation  rates  from  Eqs.  (24)  through  (31)  are 


K+  =  7.79  deg/day  (34) 

K-  =  7.11  deg/day 

<«+>  *  7.79(2  -  2.5  sin^  98’) 

*  -3.52  deg/day 

<w_>  =  7.11(2  -  2.5  sin^  98’) 

=  -3.21  deg/day 

<Q  >  *  -7.79  cos  98’ 

+ 

=  1.08  deg/day 

« 

<Q_>  =  -7.11  cos  98’ 

=  0.989  deg/day  (35) 
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The  relative  mean  apsidal  and  nodal  rotation  rates  for  the  two  groups 
of  particles  are 

<Au>  =  <u  >  -  <u  > 

+ 

=  -0.31  deg/day 

<AQ>  =  <S  >  -  <5  > 

+ 

=  0.091  deg/day  (36) 

Relative  closure  of  the  two  groups  of  particles  occurs  when  <^>  and  <a5> 
are  180*. 

are 

180“ 

S  — ■ 

I  <Aa)>  I 

180*  j  It 

*  =581  days  =  1.6  years 
180* 

1  Rn^ 

*  091= 

=5.4  years  (37) 


The  corresponding  mean  periods 


(d 


"Q 


2.10  SUMMARY  AND  CONCLUSIONS 


A  new  and  simple  method  has  been  described  in  the  report  which  can  be 
used  to  examine  the  short-term  collision  hazard  for  a  spacecraft  passing 
through  a  cloud  of  particles  resulting  from  a  breakup  of  an  object  in  orbit. 
The  method  uses  linearized  equations  for  relative  motion  to  compute  the  volume 
of  the  cloud  of  particles  as  a  function  of  the  initial  spread  velocities, 
orbital  angular  rate,  and  time. 
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A  representative  short-term  collision  hazard  for  a  spacecraft  passing 
through  the  center  of  the  cloud  was  computed  under  certain  simplifying 
assumptions.  The  results  showed  that  the  greatest  hazard  occurs  at  or  shortly 
after  the  breakup  event  when  the  cloud  volume  is  small  and  the  density  large. 
The  collision  hazard  was  found  to  decrease  rapidly  with  time. 

The  effects  of  earth's  oblateness  on  the  evolution  of  the  cloud  also 
were  examined.  The  results  showed  that  the  initially  toroidal  debris  cloud 
becomes  a  spherical  shell  due  to  nodal  and  apsidal  rotation  effects  of  the 
particle  orbits. 

The  long-term  collision  hazard  thus  is  seen  to  be  greatly  reduced  over 
a  period  of  several  years.  A  more  general  case  of  a  spacecraft  entering  an 
expanding  debris  cloud  is  considered  in  Section  8.  The  probabilities  of 
collision  with  a  particle  in  the  cloud  are  calculated  using  program  DEBRIS, 
where  the  spacecraft  path  through  the  cloud  is  determined  and  the  corresponding 
probabilities  of  collision  are  obtained. 
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3.  COMPARISONS  OF  EXACT  AND  LINEAR  SOLUTIONS 


3.1  INTRODUCTION 

Whenever  a  dynamical  problem  is  evaluated,  the  question  arises  of 
whether  to  use  the  exact  equations  of  motion  or  the  linearized  set.  The  pri¬ 
mary  concern  when  using  the  linearized  version  is  the  violation  of  the  assump¬ 
tions  made;  for  if  the  assumptions  remain  valid  over  the  range  of  interest, 
then  the  linearized  solution  should  be  an  accurate  portrayal  of  what  is  actu¬ 
ally  happening.  In  this  section,  the  dynamical  equations  of  relative  motion 
between  two  bodies  are  studied,  both  the  linearized  version  and  exact  version, 
and  the  solutions  are  compared  for  accuracy.  Also,  when  the  most  accurate  of 
these  results  is  used  to  determine  the  position  of  many  particles  relative  to 
a  rotating  reference  point,  an  approximate  volume  in  space  can  be  found  which 
could  represent  the  region  of  space  occupied  by  debris  from  an  on-orbit 
breakup  of  a  spacecraft. 

3.2  ANALYSIS 


Assume  the  satellite  in  orbit  as  a  point  moving  in  a  known  gravity 
field.  At  any  given  time,  the  position  and  velocity  vectors,  in  an  earth- 
based  inertial  frame,  are  known  from  telemetric  measurements.  When  the 
orbital  breakup  begins,  a  multitude  of  particles  will  mpve  away  from  their 
previous  position  (the  center  of  mass,  or  CM)  at  some  velocity.  This  is 
similar  to  a  probe  being  ejected  from  a  host  ship  at  a  known  initial 
velocity.  This  ejecta  now  moves  in  its  own,  unique  orbit. 


By  propagating  the  Av  vector  in  all  directions  and  in  three  dimen¬ 
sions,  one  can  determine  the  maximum  distances  from  the  former  center  of  mass. 
When  the  debris  volume  is  approximated  as  a  growing  torus,  as  shown  in 
Figure  14,  the  time-varying  volume  is  approximately  equal  to 
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where 


X  =  maximuni  distance  in  the  -x  direction  (leading  CM) 
max 

X  .  =  maximum  distance  in  the  +x  direction  (trailing  CM) 

min 

y  =  maximum  distance  in  the  +y  direction  (outside  CM) 
max 

y  .  =  maximum  distance  in  the  -y  direction  (inside  CM) 

min 

z  =  maximum  distance  in  the  +z  direction  (above  CM) 
max 

z  .  =  maximum  distance  in  the  -z  direction  (below  CM) 

min 

R  =  distance  from  the  center  of  the  earth  to  CM 


Figure  14.  Torus  Model 

Analytical  expressions,  on  the  other  hand,  allow  for  simplified 
solutions  to  complex  problems.  For  circular  orbits,  an  approximate  solution 
describing  relative  motion  between  two  bodies  revolving  about  the  same 
gravitational  attracting  mass  is  the  Clohessy-Wiltshire  equations.  In  matrix 
form,  the  analytical  solution  is 
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Since  the  particles  are  at  the  origin  of  the  coordinate  frame  at  t  =0,  the 
initial  state  vector  is 


=  {0.  0,  0,  x^, 


(^1) 


Thus,  if  position  information  is  all  that  is  desired,  then  the  system  reduces 
to 
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{X}  =  [<|>i2]  {Xo}  (43) 

where  submatrix  in  the  upper  right-hand  position  of  the  tran¬ 


sition  matrix  <I»(t,0)  and  is  similar  to  Eq.  (4).  The  volume  enclosed  in  this 
space  is  found  by  the  determinant  of  multiplied  by  4^/3  [an  expansion 
of  Eq.  (5)] 

V  =  ^  |sin  0  [(-30  +  4  sin  0)  sin  0  +  4(1  -  cos  (44) 
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where  0  =  ut,  and  Av  is  the  particle  relative  velocity. 


.  '1  ‘1  '2.  1/2 

Av  =  (x  +  y  +  z  )  (45) 

0-^0  0 

Applying  what  is  now  known  to  the  problem  of  on-orbit  breakup  of  a 
spacecraft,  one  can  determine  the  range  of  validity  of  the  Clohessy-Wiltshire 
equations . 

3 . 3  RESULTS 


The  range  of  validity  can  be  determined  by  comparing  the  effects  of  many 
different  Av  vectors  on  the  exact  and  linear  solutions  of  the  problem. 

In  each  case  studied,  the  Av  vector  was  propagated  in  six  directions:  positive 
and  negative,  tangential,  radial,  and  out  of  plane.  The  orbit  in  study  is  a 

-3 

low-altitude  orbit,  with  an  angular  rate,  «,  of  1.1  x  10  rad/sec.  Allowing 
the  magnitude  of  Av  to  range  from  1  to  1000  ft/sec  gave  a  wide  range  of  results 
which  are  summarized  below. 

A  good  comparison  between  the  exact  and  linear  solutions  is  to  look  at 
the  absolute  difference  between  the  two  solutions.  Figures  15  through  19  show 
the  differences.  The  ordinate  is  the  log  to  the  base  10  of  the  deviation  in 
nautical  miles,  while  the  abscissa  is  the  number  of  orbit  revolutions. 

Figures  15  and  16  are  deviations  in  x  versus  time  and  y  versus  time, 
respectively,  for  a  radial  ejection.  The  z  versus  time  has  been  left  out, 
since  the  out-of-plane  component  for  the  linear  solution  is  zero.  Generally, 
for  spread  velocities  less  than  100  ft/sec,  the  maximvim  deviation  is  on  the 
order  of  magnitude  of  1  nmi  after  two  orbits.  Figures  17  and  18  examine  the 
same  scenario,  except  now,  for  a  tangential  (parallel  to  the  velocity  vector) 
ejection.  For  spread  velocities  of  less  than  100  ft/sec,  the  deviation  is  on 
the  order  of  magnitude  of  10  nmi  after  two  orbits.  The  deviation  in  z  is  left 
out  for  similar  reasons  as  before. 
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Figure  19  is  the  deviation  in  z  versus  time  for  an  out-of-plane 
ejection.  For  a  spread  velocity  of  100  ft/sec,  the  deviation  grows  to  about 
0.1  nmi  in  approximately  two  orbits. 

The  comparison  of  the  torus  model  to  the  Clohessy-Wiltshire  volume 
model  is  shown  in  Figure  20.  Clohessy-Wiltshire  has  a  greater  volume, 
generally,  over  the  first  orbit.  The  torus  model  has  two  more  zero  points 
than  the  linear  model;  however,  th^  two  models  do  have  their  common  zero 
points  at  approximately  the  same  time. 

3.4  SUMMARY 

For  small  changes  in  relative  velocities  (<  100  ft/sec),  the  Clohessy- 
Wiltshire  equations  prove  to  be  adequate  for  at  least  the  first  couple  of 
orbit  revolutions.  However,  as  time  increases,  the  equations  degrade  as  the 
nonlinear  terms,  previously  neglected,  begin  to  have  a  larger  effect.  When 
large  velocity  changes  are  considered  (>  1000  ft/sec),  the  radial  and 
tangential  position  components  for  orbital  plane  ejections  begin  to  deteriorate 
quickly  as  time  progresses.  However,  the  out-of-plane  component  of  position 
for  an  out-of-plane  ejection  was  quite  accurate  for  both  the  linear  and  exact 
equations . 

When  considering  large  time  periods  (days),  however,  one  must 
necessarily  use  the  exact  equations,  since  the  linear  equations  are  good  only 
for  a  short  time  following  the  ejection. 
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Figure  15.  Deviation  Between  Exact  and  Linear  Solution: 
X  versus  Time;  Radial  Ejection 
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Figure  16.  Deviation  Between  Exact  and  Linear  Solution 
y  versus  Time;  Radial  Ejection 


DEVIATION  (nmi)  ^  DEVIATION  (nmi 


NUMBER  OF  ORBITS 


figure  17.  Deviation  Between  Exact  and  Linear  Solution: 
X  versus  Time;  Tangential  Ejection 
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Figure  18.  Deviation  Between  Exact  and  Linear  Solution: 
y  versus  Time;  Tangential  Ejection 
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VOLUME  DIMENSIONLESS  DEVIATION  (nmi) 
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Figure  19.  Deviation  Between  Exact  and  Linear  Solution: 
z  versus  Time;  Out-of-Plane  Ejection 


Figure  20.  Volume  versus  Timet  Clohessy-Wiltshire 
and  Torus  Approximation 
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4.  THE  EFFECTS  OF  PERTURBATIVE  FORCES 
ON  THE  DEBRIS  CLOUD  EVOLUTION 


4.1  INTRODUCTION 

In  previous  sections,  the  volume  of  a  debris  cloud  following  an  on-orbit 
breakup  of  a  spacecraft  assiuned  an  inverse-sqxiare  gravitational  field  and 
failed  to  include  any  additional  external  forces  acting  upon  the  individual 
debris  particles.  In  actuality,  there  are  always  external  forces  acting  in 
addition  to  higher  order  harmonic  terms  in  the  gravity  field.  This  section 
reexamines  the  effects  of  one  such  harmonic  term,  (earth  oblateness 
term),  and  a  potentially  important  external  force,  atmospheric  drag,  on  the 
shape  and  size  of  a  debris  cloud.  The  approach  used  is  based  upon  a 
simplified,  first-order  formulation  which  provides  a  cursory  but  timely 
estimate  of  the  propagation  time  for  the  cloud  to  cover  the  entire  globe  as 
well  as  an  early  estimate  of  volume  effects  used  for  density  calculations. 

4.2  BACKGROUND 


As  shown  in  Sections  2  and  3,  the  linearized  Clohessy-Wiltshire 
equations  provide  an  adeqxiate  solution  for  the  relative  motion  between  two 
objects  in  space  with  low  initial  relative  velocities.  However,  there  exist 
solutions  to  the  linearized  volume  equations  that  produce  a  "pinch  point"  or  a 
point  of  zero  volume.  This  point  occurs  at  integer  number  of  revolutions  after 
the  orbital  breakup  begins.  Because  of  this  pinch  point,  the  density  of  debris 
becomes  infinite  (number  of  particles  divided  by  the  volume)  and  thus  produces 
an  unacceptable  collision  hazard.  However,  since  this  volume  is  initially 
infinitesimal,  there  is  very  little  chance  of  a  collision  with  another 
spacecraft. 

In  the  true  situation,  there  will  not  be  an  exact  regrouping  of  the 
debris.  External  forces  plus  variations  in  the  earth's  gravity  field  will 
cause  the  volume  to  be  perturbed,  thus  relieving  the  singularity  in  density  at 
integer  number  of  revolutions.  Therefore,  accounting  for  some  of  these 
perturbations  will  provide  a  better  understanding  of  what  actually  happens  to 
this  debris  cloud. 
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4.3 


ANALYSIS 


As  shown  in  Section  2,  the  dimensionless  volume  of  a  debris  cloud  is 
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where  9  is  the  angle  swept  out  by  the  rotating  reference  frame  with  respect 
to  the  inertial  frame.  The  matrix  [M]  is  the  state  transition  matrix  relating 
relative  position  to  initial  velocity. 


In  this  section,  the  elements  of  [M]  are  slightly  perturbed  by  and 
atmospheric  drag,  and  the  results  are  compared  to  those  obtained  in  the 
linearized  solution. 


The  first  step  is  to  perturb  the  functions  a 
formulation  is 


One  possible 


a^j^  =  -39  +  4  sin  9 

3^2  “  ~®21  “  ^  ~  f^(t)] 

a22  =  sin  9[1  -  f2(t)  +  f2(t)/|sin  9|] 

^33  “ 
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These  functions  have  initial  values  of  zero  and  steady-state  values  of 


f,(tf)  =  1 


f2(tp  =  1 


^  a  sin  1 


where  i  is  the  inclination  of  the  target  spacecraft,  u  is  the  rate  of 
increase  of  9,  and  t^  is  the  ending  time  of  the  approximation. 


Assume  a  nearly  linear  relationship  for  the  functions  fj^(t),  f^Ct), 


and  g(t),  such  that 


f^(t)  =  C^t  +  C^(9  -  sin  6) 
f2(t)  =  C2t  +  Cj(9  -  sin  0) 

g(t)  =  C3t  (50) 

where  the  constants  i=l,2,...,5  are  to  be  determined,  and  depend  upon 
the  orbital  elements  (altitudes  at  apogee  and  perigee,  and  inclination)  of  the 
debris  source. 

The  functions  fj^(t)  and  f2(t)  spread  the  debris  out  in  the  radial 
direction;  they  come  directly  from  the  spreading  due  to  atmospheric  drag  and 
the  rotation  of  the  line  of  apsides  of  the  individual  particle  orbits.  The 
function  g(t)  is  an  out-of-orbit  plane  change,  this  time  due  to  a  shifting  in 
the  right  ascension  of  ascending  nodes  of  the  debris  particle  orbits. 

However,  the  right  ascension  of  ascending  nodes  rotates  about  the  polar  axis 
of  the  earth,  and  not  the  normal  axis  to  the  orbit. 


47 


4.4  EFFECTS 


The  constants  C^,  C^,  and  are  due  to  the  imperfections  in  the 
earth's  gravitational  field;  i.e.,  the  zonal  harmonic  term  J^.  The  and 
terms  relate  the  shifting  in  the  line  of  apsides,  while  the  term 
relates  the  shift  in  the  line  of  nodes  of  the  orbit. 

First,  take  the  equation  for  the  shifting  of  the  line  of  apsides  at  an 
altitude  h  as 


CO  = 


9.96399 


(1  - 


2.2 

e  ) 


R 


(R  +  h) 


3-5 


[2  -  ^  sin^(i)] 


deg 


mean  solar  day 


(51) 


This  is  with  respect  to  an  earth-centered,  inertial  reference  frame.  Assume 
as  a  center  of  a  rotating  reference  frame  (with  respect  to  the  inertial  frame) 
the  position  of  an  object  before  it  breaks  up  in  orbit.  After  the  breakup 
occurs,  groups  of  particles  will  be  traveling  with  a  greater  co  than  the 
satellite,  had  the  breakup  not  occurred.  Likewise,  there  are  groups  of 

particles  with  smaller  co  values.  Statistically,  by  finding  the  average 

•  * 

greater  co  and  the  average  lesser  co,  one  can  find  a  relative  shifting  in 
the  line  of  apsides.  For  a  greater  co,  the  semi-major  axis  decreases.  This 
probabilistic  change  (Aa)  in  the  semi-major  axis,  a,  is  a  function  of  the 
change  in  velocity  as  per  Eq.  (28) 

<Aa>  =  a  (Av/v)  (52) 

where  v  is  the  orbital  velocity  at  breakup  with  respect  to  an  earth-centered 
inertial  reference  frame.  The  relative  shift  in  the  line  of  apsides  for  a 
circular  orbit  is 

<Aco>  =  AK  [2  -  I  sin^(i)]  (53) 

where 

AK  =  K+  -  K_ 
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and,  converting  to  appropriate  units 
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and 
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Similarly,  for  the  rotation  of  the  line  of  nodes,  the  equation 


-9.9639  /R 

/I  2.2U/ 
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deg 


mean  solar  day 


(56) 


relates  the  change  in  nodes  as  a  function  of  orbital  elements. 


Again,  we  are  interested  in  a  relative  rotation  in  the  line  of  nodes,  so 


<a5>  =  I AK  cos(i) I 


(57) 


The  time  that  it  takes  the  line  of  nodes  to  shift  180®  is 


and  the  time  it  takes  the  line  of  apsides  to  shift  180®  is 

T  - 

<A(i)> 


(58) 


(59) 
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The  constants  C^,  C2»  and  C3  are  therefore 


=  r 

(1) 


CO 


=3 ' 1;^ 


(60) 


where  sin(i)  is  called  a  volume  limit  factor  with  limits  on  the  functions  as 

fl(t)  =  Cit  <  1 
f2(t)  =  C2t  <  1 


(1) 


(61) 


4.5  ATMOSPHERIC  DRAG 

The  effects  of  atmospheric  drag  affect  only  the  radial  and  tangential 
components  of  the  position.  This  is  an  additive  effect  and  is  unrelated  to 
J^.  We  can  now  assume  a  new  expression  for  the  perturbation  functions 
fj^(t)  and  f2(t).  The  function  g(t)  remains  unchanged,  since  this  is  the 
out-of-plane  function.  As  in  Eq.  (50),  assume 

f,(t)  =  C, t  +  C,  (0  -  sin  9) 

1  14 


f2(t) 


f^(t) 


The  constant  C(^  is  now  to  be  determined. 


(62) 
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From  Reference  7,  the  drop  in  altitude  due  to  atmospheric  drag  is 


Ah 


_Jie£_ 

W/C„A 


(0  -  sin  0) 


(63) 


where 

P 

g 

W 

A 


gravitational  constant 
atmospheric  density 

acceleration  due  to  gravity  at  sea  level 
weight  of  debris  particle 
drag  coefficient 

cross-sectional  area  of  debris  particle 


For  a  group  of  particles  with  a  certain  Av,  assuming  an  isotropic 
expansion,  half  of  the  particles  would  move  into  a  higher  energy  orbit,  and 
half  would  move  into  a  lower  energy  orbit.  Therefore,  the  immediate  conse¬ 
quence  is  that  the  particles  with  a  lower  energy  orbit  will  be  affected  more 
by  the  atmosphere  than  those  at  a  higher  energy  level.  Multiplying  by  a 
statistical  factor  of  0.5  and  nondimensional zin?  •  Eq.  (63)  becomes 


Ah' 


(0  -  sin  0) 


(64) 


where  Ah'  is  a  nondimensional  drop  in  altitude.  The  factor  C4  is  the 
average  over  one  revolution 
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where  t=  VI/C^A  (ballistic  coefficient). 


A  useful  expression  for  the  atmospheric  density  is  (Ref.  8) 


p  =  Ke 


-h/c 


(66) 
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where 


K  =  reference  density  =  4. 80711  x  10  ^  Ib/ft^ 
h  =  altitude 

c  =  reference  altitude  =  12.58  nmi 

The  altitude  function  h  is  dependent  upon  9,  while  everything  else  is 
approximately  independent  of  9.  By  definition 

h  =  r  -  R  =  a.,(l  -  e„  cos  9)  -  R  (67) 

N  N 

The  values  a^^  and  are  semi-major  axis  and  eccentricity  of  the  new  lower 
energy  orbit,  respectively,  with 


(68) 

(69) 

where  a  is  the  circular  radius  of  the  object  orbit  that  is  breaking  up. 
Assembling  the  various  expressions,  Eq.  (65)  becomes 


a^j  =  a  [1  -  (Av/v)] 

and 

^  _  (Av/v) 

~  1  -  (Av/v) 


2ir 


C  -  1 _ Q-5pg 

-  2it  /Av  .  g 

c 


(^) 


K  f  exp{-[ajj(l  -  ej^  cos  9)  -  R]/c}d9 


(70) 


Integrating  Eq.  (70) 


‘ih 


3/2 


O.SpKg 

c 

\  CJ/  c 

a(Av/v) 

1/2 


exp{[R  -  a^(l  -  ej^)]/c} 


(71) 


and,  as  before. 


c 

^5  ■  2 
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4.6 


RESULTS 


Take,  for  example,  the  following  two  cases. 

Case  1:  Assume  that  the  object  breaking  up  is  at  an  altitude  of 

200  nmi,  inclined  at  45*  with  respect  to  the  equator.  Figures  21  and  22  show 

the  constants  C^,  i  =  1,2.. .5,  for  an  applied  Av.  Assuming  a  Av  of  100  m/s, 

we  note  that  Figure  23  shows  the  volume  versus  time  profile  with  only,  drag 

only,  no  perturbations,  and  drag  and  J-.  At  9  =  2Tr,  the  volume  (in  non- 

-3  ^  3 

dimensional  form)  is  4.7  x  10  (diraensionalizing,  V  =  3200  km  ) .  In  this  case, 
it  will  take  322  days  for  the  line  of  apsides  to  shift  180°,  and  5.7  years  for 
the  cloud  to  envelop  the  earth. 

Case  2;  Again,  we  have  a  circular  orbit  with  an  altitude  of  500  nmi, 

with  an  inclination  of  45°.  Figures  24  and  25  show  the  constants  C^,  i  =  1, 

2... 5  for  breakup.  Figures  26  show  the  volume  versus  time  profile  with 

only,  drag  only,  no  perturbations,  and  J.  drag.  At  0  =  2ir,  the  nondimensional 
-5  ^3 

volume  is  6.5  x  10  (dimensionally,  V  =  63  km  )  which  produces  a  finite  density 
at  0  =  2-rr.  For  this  example,  it  will  take  410  days  for  the  line  of  apsides  to 
rotate  180°,  and  7.2  years  for  the  cloud  to  envelop  the  earth. 

4.7  SUMMARY 


The  addition  of  atmospheric  drag  into  the  volumetric  equations  spreads 
the  cloud's  shape  into  a  plane  at  the  integer  orbit  pinch  points.  Only  with 
the  addition  of  do  we  have  a  spreading  out  in  the  out-of-orbit  direction, 
thus  creating  a  volume  at  the  pinch  points.  Although  the  volumes  at  the  pinch 
points  in  the  first  few  revolutions  are  small,  they  are  not  zero  and  thus  yield 
a  more  accurate  assessment  of  the  potential  collision  hazard.  Also,  as 
expected,  for  a  higher  altitude  orbital  breakup,  the  atmospheric  drag  has  less 
of  an  effect  than  for  a  lower  orbit  breakup. 
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DRAG  PARAMETERS  (non-dimensional) 


Figure  21.  J2  Parameters,  C^,  C2»  C3,  versus  Av 
200-nmi  Circular  Orbit 


Figure  22.  Drag  Parameters,  C4,  C5,  versus  Av; 
200-nmi  Circular  Orbit 


5 


NUMDER  ()F  ORIIITS 


Figure  23a.  Volume  versus  Time  for  Breakup  at  200-nmi  Altitude 
J2  Only 


NUMBER  OF  ORBITS 


Figure  23b.  Volume  versus  Time  for  Breakup  at  200-nmi  Altitude 
Drag  Only 


Figure  23c.  Volume  versus  Time  for  Breakup  at  200-nmi  Altitude, 
No  Perturbations 


Figure  23d.  Volume  versus  Time  for  Breakup  at  200-nmi  Altitude, 
J2  and  Drag 
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Figure  2k 
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J2  Parameters,  Cj,  C2,  C3,  versus  Av 
500-nmi  Circular  Orbit 
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Figure  25.  Drag  Parameters,  C4,  C5,  versus  Av 
500-nmi  Circular  Orbit 
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Figure  26a.  Volume  versus  Time  for  Breakup  at  500-nmi  Altitude 
J2  Only 
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Figure  26b.  Volume  versus  Time  for  Breakup  at  500-nmi  Altitude 
Drag  Only 
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Figure  26c.  Volume  versus  Time  for  Breakup  at  500-nmi  Altitude 
No  Perturbations 


NUMBER  OF  ORBITS 


Figure  26d.  Volume  versus  Time  for  Breakup  at  500-nmi  Altitude 
J2  and  Drag 


5.  COMPUTING  FRAGMENT  DIRECTIONS  AND  ORBITS  AFTER  COLLISION 


5.1  INTRODUCTION 


Section  6  describes  a  method  for  determining  the  masses  and  velocities 
of  fragments  that  result  from  a  hypervelocity  collision.  The  velocities 
derived  are  relative  to  the  center  of  mass  (CM).  This  section  calculates  the 
center  of  mass  velocity  vector  and  determines  the  fragments'  directions 
relative  to  the  center  of  mass.  Both  sets  of  results  are  then  used  to  deter¬ 
mine  the  fragments'  orbits.  An  example  is  used  to  illustrate  the  method.  The 
results  are  tabularized  and  graphed. 

5.2  METHOD  OF  ANALYSIS  AND  EQUATIONS 

There  are  three  steps  in  the  methods  of  analysis  to  determine  the 
fragment  orbits:  (a)  to  determine  the  orbit  of  the  center  of  mass  of  the  two 
objects  colliding;  (b)  to  calculate  the  masses,  velocities  (Section  6),  and 
directions  of  the  fragments  relative  to  the  center  of  mass;  and  (c)  to  input 
the  above  results  into  an  orbit  element  conversion  program  to  compute  orbital 
properties  of  the  fragments.  Figure  27  summarizes  these  three  points  and 
illustrates  the  center  of  mass  orbit  of  two  intersecting  orbits. 

5.2.1  Coordinates  of  Center  of  Mass 


To  find  the  coordinates  for  the  center  of  mass  of  the  two  objects,  let 
be  the  position  vector  of  the  two  objects  at  intersection,  and  let  Vj^  and  v^  be 
their  respective  velocity  vectors  (as  shown  in  Fig.  27)  such  that 


=  ri  +  rj  +  rk 

(72) 

X  y  z 

1  =  ''lx‘  *  ''lyj  *  ''u'' 

(73) 

2  -  ''2x‘  *  >5  ♦  ’2.^ 

(7^.) 

The  conservation  of  angular  momentum  is  then  applied  in  which  the  total 
angular  momentum  of  the  system  at  the  intersection  point,  is  equal  to  the 
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METHOD  OF  ANALYSIS 


•  DETERMINE  CENTER  OF  MASS  ORBIT  OF  TWO  VEHICLES 

•  CALCULATE  VELOCITIES  AND  DIRECTIONS  OF  FRAGMENTS 
RELATIVE  TO  CENTER  OF  MASS 

•  DETERMINE  ORBITAL  ELEMENTS  OF  DEBRIS 


r  =  radius  vector  of  encounter  point  V2  =  velocity  vector  of  Object  2 

CM  =  center  of  mass  of  Objects  1  and  2  Vqi^  =  velocity  vector  of  center  of  mass 
Vi  =  velocity  vector  of  Object  1 


Figure  27.  Center  of  Mass  Vector  from  Two  Intersecting  Objects 
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sum  of  the  angular  momenta  of  the  two  objects,  ^2'  quantity,  P,  refers 

to  the  total  angular  momentum  as  if  it  were  concentrated  at  the  center  of  mass. 
This  is  precisely  what  we  want,  since  we  are  trying  to  find  the  velocity  of  the 
center  of  mass.  Using  the  conservation  principle,  we  get 


I  . 


(75) 


where 


=  r  X 


P2  =  r  X  m^v^ 


P  =  r  X  Mv, 


CM 


M  =s  m^^  +  m^ 


v_„  =vi+vj+vk 
CM  X  y'  z 


(center  of  mass  velocity) 


(76) 

(77) 

(78) 

(79) 

(80) 


Solving  Eq,  (75),  using  the  rest  of  the  relations,  yields  the  components 
of  the  center  of  mass  velocity  vector 


Vx  =  (raivix  +  in2V2x)/M 


(81) 


Vy  =  (m^viy  +  ro2V2y)/M 


(82) 


(83) 


Thus,  the  center  of  mass  velocity  is  expressed  in  terms  of  mass  of  the 
two  objects  and  their  velocities. 


5.2.2  Fragment  Velocity  Relative  to  the  Center  of  Mass 

The  kinetic  energy  of  the  colliding  objects  and  the  center  of  mass  are 
used  to  compute  the  fragment  velocities. 
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The  total  kinetic  energy  of  the  two  objects  before  collision  is  given  by 


2  2 

KE  (before  collision)  = 


After  the  collision,  the  total  kinetic  energy  of  the  system  can  be 
expressed  by  the  sum  of  the  kinetic  energy  of  the  total  mass  moving  with  the 
velocity  of  the  center  of  mass  and  that  due  to  the  motions  of  the  individual 
particles  relative  to  the  center  of  mass  (Ref.  9).  The  following  equation 
expresses  this  as 


KE  (after  collision)  = 


where 


mj[  =  mass  of  i-th  fragment 

p  =  velocity  magnitude  of  i-th  fragment 
relative  to  the  center  of  mass 

n  =  number  of  fragments 

If  the  collision  is  assumed  to  not  conserve  kinetic  energy,  then 
Eqs.  (84)  and  (85)  can  be  related  in  Eq.  (86)  (canceling  the  1/2's) 


1  =  1 

where  KE^^^g^  =  twice  the  amount  of  kinetic  energy  lost  in  the  collision. 


The  unknowns  in  Eq.  (86)  are  m^,  p^,  and  n.  Thus,  rearranging  Eq.  (86) 


yields 


=  m^v^2  +  m^v^^  -  (m^  +  m2)v2^  -  KE^ 
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As  mentioned  earlier.  Section  6  describes  a  method  for  determining  the 
number  of  fragments,  n,  the  fragment  masses,  m^,  and  velocities,  p^,  that 
solve  Eq.  (87). 

5.2.3  Determining  Fragment  Directions 

Conservation  of  momentum  is  used  to  determine  fragment  direction.  A 
random  direction  relative  to  the  center  of  mass  was  given  to  the  first  frag¬ 
ment.  The  second  fragment  was  given  a  direction  (relative  to  the  center  of 
mass)  opposite  the  first.  In  this  way,  the  momentums  imparted  to  the  first 
two  fragments  relative  to  the  center  of  mass  cancel.  This  assumes  that  these 
fragments  have  the  same  mass  and  velocity.  The  process  is  repeated  until  all 
fragments  have  been  assigned  directions.  The  directions  are  chosen  randomly 
from  the  uniform  directions  derived  in  Appendix  B. 

5.2.U  Determining  Fragment  Orbits 

After  completion  of  the  previous  analysis  in  Sections  5.2  and  6,  the 
center  of  mass  position  and  velocity  should  be  known  and  the  fragments  should 
all  have  assigned  velocities  and  directions  relative  to  the  center  of  mass. 

The  velocity  of  the  center  of  mass  and  the  fragment  velocities  relative  to  it 
can  be  vectorial ly  added  to  obtain  the  inertial  velocity  of  each  of  the 
fragments.  At  this  point,  orbit  element  conversions  are  used  to  convert  the 
elements  from  earth  centered  inertia  (ECI)  to  other  coordinate  systems  in 
order  to  determine  the  orbital  properties  of  the  fragments.  In  this  analysis, 
the  On-Line  Orbital  Mechanics  (OLOM)  program  (Ref.  10)  was  used  to  perform  the 
conversions . 

5.3  DATA  GENERATION 


Once  the  orbit  elements  of  the  fragments  have  been  determined,  the  data 
can  be  presented  in  various  ways.  In  our  analysis  we  presented  data,  typical¬ 
ly,  in  three  formats.  The  first  is  a  cable  showing  the  percentage  of  fragment 
perigees  above  100  nmi,  the  percentage  between  0  and  100  nmi,  and  the  percen¬ 
tage  below  the  earth's  surface.  Thus,  the  table  shows  how  many  fragments  will 
remain  in  orbit  and  reenter  sometime  later,  and  those  that  will  reenter 
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immediately.  An  example  is  shown  in  Table  3.  In  this  example,  the  colliding 
bodies  have  a  mass  ratio  of  15.  The  heavier  object  is  in  low  earth  orbit,  and 
the  relative  velocity  at  collision  is  23,183  ft/sec.  The  first  three  columns 
are  the  result  of  the  analysis  in  Section  6. 

The  second  way  to  present  data  is  to  plot  the  fragment  apogee  and  peri¬ 
gee  altitudes  versus  the  fragment  periods.  An  example  of  this  plot  is  shown 
in  Figure  28.  The  0  nmi  altitude  represents  the  earth's  surface.  The  figure 
resembles  two  wings  that  meet  at  a  point.  The  upper  wing  shows  the  apogee 
altitudes,  and  the  lower  wing  shows  the  perigee  altitudes  versus  the  periods. 
Thus,  two  points  are  plotted  for  each  fragment.  Note  the  correspondence 
between  the  perigees  plotted  in  Figure  28  and  the  percentages  listed  in 
Table  3.  The  APL  graphics  package,  EZPLOT  (Ref.  11),  was  used  in  conjunction 
with  OLOM  to  generate  Figure  28.  It  would  be  possible  also  to  plot  any  of  the 
fragment  orbital  elements. 

Notice  in  Figure  28  that  the  right  side  of  the  upper  wing  and  the  left 
side  of  the  lower  wing  are  not  straight.  They  have  a  slight  bow  in  them. 

This  bow  is  only  noticeable  with  the  high  spread  velocities.  This  example 
used  spread  velocities  of  over  5000  ft/sec.  In  Figure  B-4  of  Appendix  B, 
these  parts  of  the  wings  appear  straight;  the  spread  velocity  used  for 
Appendix  B's  example  is  1000  ft/sec.  If  the  spread  velocity  is  large  enough, 
the  bow  appears  regardless  of  whether  the  fragment  spread  is  from  a  circular 
or  an  eccentric  orbit  and  regardless  of  the  position  within  the  eccentric 
orbit . 


The  third  method  of  presenting  data  is  to  plot  footprints  that  show  the 
region  on  the  earth  over  which  reentering  debris  will  fall.  The  ECI  coordi¬ 
nates  of  the  reentering  fragments  are  input  to  the  program  PECOS  (Parametric 
Examination  of  the  Cost  of  Orbit  Sustenance,  Ref.  12)  and  propagated  until  a 
specified  altitude  is  reached.  PECOS  will  print  out  the  latitude  and  longitude 
of  the  fragments  at  the  specified  altitude.  Any  program  can  be  used  that  can 
accurately  propagate  the  orbital  elements  under  the  influence  of  drag.  The 
latitudes  and  longitudes  are  then  plotted  onto  a  map  using  EZPLOT. 
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Mass  ratio  is  of  colliding  bodies 
Av  is  relative  to  the  center  of  mass 
Spacecraft  orbit:  Low  altitude  circular 
Relative  velocity  of  collision:  23183  ft/sec 


Figure  28.  Apogee  and  Perigee  Altitudes  of  Each  Fragment 
versus  Its  Period;  MR  =  15;  NF  =  992 


Since  the  concern  of  our  analysis  was  the  extent  of  the  footprint  and 
not  the  distribution  of  fragments  within  it,  it  was  not  necessary  to  propagate 
the  elements  of  all  the  fragments.  In  this  analysis,  only  seven  fragment 
directions  (per  mass  group)  that  would  define  the  contour  of  the  footprint 
(for  that  particular  group)  are  propagated.  All  the  directions  used  lie  in 
the  plane  defined  by  the  velocity  vector  of  the  center  of  mass,  v„„,  and  the 
cross  product  of  the  position  and  velocity  vectors,  rxv^^^  (Fig.  29).  More 
specifically,  the  directions  chosen  were  in  the  uprange  half  of  the  figure  (180 
to  360® — the  shaded  region). 
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Figure  30  shows  a  typical  footprint  that  includes  three  contours. 

There  is  one  contour  for  each  mass  group  (0.1,  1,  and  10  lb).  Each  contour 
was  formed  by  connecting  seven  "x's."  The  seven  x's  are  the  latitude/ 
longitude  "impact"  locations  that  result  from  the  seven  fragment  directions 
propagated  by  PECOS,  The  three  contours  are  centered  over  the  groundtrack  of 
the  orbiting  object.  The  10-lb  and  heavier  fragments  would  be  contained 
within  the  smallest  contour,  1.0-lb  and  heavier  in  the  middle  contour,  and 
0.1-lb  and  heavier  in  the  largest  contour. 

Note  that  the  figure  does  not  show  how  many  of  the  fragments  will 
survive  reentry  (if  any)  and  that  the  masses  listed  are  the  masses  before 
reentry,  not  at  impact.  Thus,  even  though  the  curves  show  regions  over  which 
fragments  of  various  "preentry"  masses  may  impact,  no  fragments  may  actually 
reach  the  ground.  Further  analysis  in  this  area  is  needed  in  order  to  assess 
the  actual  threat  on  the  ground  to  impacting  debris. 

5.4  SUMMARY 


This  section  explained  the  three  steps  used  to  determine  the  orbital 
elements  of  the  fragments  resulting  from  two  colliding  objects.  First, 
calculate  the  center  of  mass  velocity  vector  from  equations  presented.  Then 
determine  the  fragment  masses,  velocities,  and  directions  relative  to  the 
center  of  mass.  Finally,  use  the  results  from  the  first  two  steps  to 
determine  the  orbit  properties  of  the  fragments. 

Three  methods  of  presenting  the  data  also  were  shown:  (a)  tabular  data 
showing  the  percentages  of  fragments  that  reenter  or  stay  in  orbit,  (b)  graphs 
of  the  orbital  elements  that  show  the  distribution,  and  (c)  a  footprint 
showing  the  regions  on  the  earth  in  which  fragments  would  be  expected  to 
impact  (if  they  survive  reentry). 
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VCM 

90°  (downrange) 


Figure  29.  Plane  Defined  by  and  t  x  Vectors  (directions 
typically  used  in  PECOS  propagation  to  generate 
fragment  footprint  lie  in  shaded  region) 
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Figure  30.  Reentering  Debris  Footprint  for  Colliding 
Objects  with  Mass  Ratio  of  15 


70 


6.  DEBRIS  GENERATION  AT  HYPERVELOCITY  COLLISION  IN  SPACE 


6.1  INTRODUCTION 


Increase  in  debris  as  a  result  of  collisions  and  explosions  in  space  has 
been  recorded.  Some  of  the  collisions  and  explosions  are  involuntary  and  some 
are  deliberate.  The  debris  can  cause  damage  to  other  satellites  in  space. 

The  number  and  velocity  of  fragments  generated  by  collision  between  two 
objects  in  space  depends  on  the  relative  velocity  of  impact;  the  masses  of 
colliding  objects;  and  the  material  properties  such  as  density,  elastic  limit, 
melting  point,  ultimate  stress  coefficient,  and  the  impact  scenario.  For 
example,  in  the  case  of  hypervelocity,  catastrophic  impact,  both  objects  will 
fragment.  In  some  other  cases,  there  will  be  a  crater  formed  in  one  object, 
whereby  the  ejected  mass  will  fragment  leaving  the  rest  of  the  object  intact. 

In  the  case  of  oblique  impact,  a  part  of  the  object  will  shear,  leaving  the 
remaining  portion  intact.  In  this  section,  only  the  hypervelocity, 
catastrophic  impact  will  be  discussed. 

This  section  presents  a  method  for  determining  the  number  of  fragments, 
their  masses,  and  spread  velocities  created  as  a  result  of  collision  between 
two  objects  in  space. 

An  analysis  of  debris  generation  at  hypervelocity  collision  in  space  has 
been  performed  assuming  conservation  of  momentum  and  energy.  A  methodology  has 
been  developed  to  determine  the  mass  of  fragments  and  their  spread  velocities 
about  the  center  of  mass.  Size  of  a  fragment  with  a  given  mass  has  been 
calculated  assuming  certain  shape  and  density. 

6.2  ANALYSIS 


When  two  masses  collide  at  near-hypervelocity  or  hypervelocity,  internal 
pressurization  develops  which,  in  turn,  generates  internal  shock  waves. 
Nebolsine,  et  al.  (Ref.  13)  have  studied  the  kinetic  energy  mechanism  at 
Physical  Sciences,  Incorporated  (PSI).  Figure  31  shows  the  shock  wave 
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generation  at  impact.  The  expressions  for  pressure,  velocity,  and  impulse  of 
the  shock  wave  are  given  below: 
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where  P  is  pressure,  v  is  the  velocity  of  shock  wave,  I  is  impulse  per  unit 
s  s  s 

area,  p  is  density  of  material,  v^  is  relative  velocity  of  impact,  d  is  dia¬ 
meter  of  projectile,  and  R  is  radius  of  shock  wave  front  which  increases  with 
time.  Shock  wave-related  parameters  are  shown  in  Figure  31.  For  fragmenta¬ 
tion,  I  >  o  t/a^  where  a^  is  the  speed  of  sound,  t  is  the  thickness  of 
s  u  t  t 

material,  and  o  is  the  ultimate  stress  coefficient, 
u 


TARGET 


Figure  31.  Shock  Wave  Propagation 
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Consider  a  case  where  two  masses  and  tn^  moving  at  velocities  v^  and 

v^  collide.  From  the  conservation  of  momentum  we  get 

m  V  +  m  V  =  (ra  +  m  )v^ _  (89) 

11  /  2  1  2  cm 


where  v  is  the  velocity  of  the  center  of  mass, 
cm 


Conservation  of  kinetic  energy  will  yield  the  following  equation: 

1  +  1  =  l(»,j  .  *  Q  (90) 

where  Q  is  the  available  energy  for  the  spread  of  fragments  created  by  impact. 
The  term  Q  is  related  to  fragment  mass  and  spread  velocity  in  the  following 
manner: 


Q  =  I  Im  ^vj  +  E'  (91) 

i 

where  E'  is  the  loss  in  energy  due  to  heat  and  light  generated  at  impact,  v^ 
is  the  spread  velocity  of  i-th  fragment  with  respect  to  the  center  of  mass, 
and  m^  is  the  mass  of  that  fragment.  The  center  of  mass  is  moving  with 
velocity  at  the  same  time,  fragments  are  spreading  with  respect  to  the 

center  of  mass.  For  the  case  v'hen  E'  =  0,  Q  (the  energy  available  to  spread 
the  fragments  about  the  center  of  mass)  depends  on  the  relative  velocity  of 
impact  and  the  ratio  between  the  two  colliding  masses  and  is  given  by 

Q  =  F  •  (KE)^^^  (92) 

2 

where  (KE)  ,  =  1/2  m,v  ,,  v  ,  =  relative  velocity  at  the  time  of  collision, 
rel  1  rel  rel  ^ 

Fraction  F  is  the  ratio  between  Q  and  the  relative  kinetic  energy,  (KE)^^j^, 
and  is  related  to  mass  ratio  R(R  =  as  shown  in  Figure  32  (Ref.  13). 

For  a  given  mass  ratio,  one  can  obtain  F  and  then  use  Eq.  (92)  to  determine  Q. 
Once  the  value  of  Q  is  determined,  the  problem  is  to  determine  fragment  mass 
and  assign  the  spread  velocity  to  satisfy  Eq.  (91).  This  is  discussed  in  the 
following  sections. 
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Figure  32.  Inelastic  Collison  Energy  Fraction  F 
versus  Mass  Ratio  R 

6.2.1  Determination  of  Fragment  Number 

To  determine  the  number  of  fragments  of  a  given  mass,  the  equation 
obtained  by  Kessler,  et  al.  (Ref.  1)  has  been  used.  This  equation  is 
empirical  and  was  obtained  on  the  basis  of  fragments  generated  in  some  test 
experiments.  This  equation  is  given  by 

N  =  k(S-)"  (93) 

where  N  is  the  number  of  fragments  of  mass  m  and  greater;  M^  is  the  total  mass 
fragmented  after  impact,  n  =  -0.75;  and  K  is  a  constant  which  depends  on  the 
rigidity  of  the  material.  For  a  rigid  material,  K  =  0.4;  for  a  nonrigid 
material,  K  =  0.8.  A  rigid  material  will  fragment  in  pieces  with  a  larger  mass, 
and  a  nonrigid  material  will  fragment  in  pieces  with  a  smaller  mass  under  the 
same  impulse.  For  collision  between  satellites,  NASA  has  determined  that  K  = 
0.45  is  a  good  value  to  be  used,  because  there  is  a  fairly  good  agreement 
between  predicted  and  observed  results  if  this  value  of  K  is  used.  Equation 
(93)  can  be  used  to  obtain  the  number  of  fragments  and  their  masses  generated 
by  hypervelocity  impact  and  not  by  explosion. 
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Figure  33  has  been  taken  from  Reference  1  and  is  based  on  the  result  of 
a  laboratory  experiment.  In  this  experiment,  a  237-gm  projectile  was  used  to 
impact  a  target  of  25  kg  mass  at  a  relative  velocity  of  3.3  km/sec.  The 
cumulative  number  of  fragments  versus  the  diameter  of  the  fragment  is  shown  in 
this  plot. 

The  fragments  generated  by  explosion  have  a  different  distribution  as  a 
function  of  size  than  the  fragments  generated  by  collision.  The  following 
equations  can  be  used  to  determine  the  number  of  fragments  generated  by 
explosion  or  breakup*: 

-U  1/2 

N  =  1.71  X  10  M^exp(-0. 02056  m  )  for  m  >  1936  gm  (94a) 

and 

-U  1/2 

N  =  8.69  X  10  M^exp(-0. 05756  m  )  for  m  <  1936  gm  (94b) 

where  is  the  mass  of  target  in  grams,  M  s  the  mass  of  fragment  in  grams, 
and  N  is  the  cumulative  number  of  fragments  with  mass  greater  than  M  grams. 


10"^  10"’  1  10  10^ 


DIAMETER  (cm) 

Figure  33.  Number  of  Fragments  Produced  from 
237-gm  Projectile 

*See  Ref.  18. 
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size  of  the  fragment  will  depend  on  the  mass  density.  For  satellite 

3 

structures,  the  density  could  vary  between  0.1  to  5  gm/cm  .  It  can  be 
assumed  that  the  small-sized  fragments  (0.1  cm  diameter  and  smaller)  will  be 
spherical  in  shape;  medium-sized  fragments  (between  0.1  and  20  cm  diameter) 
will  be  cubical  in  shape;  and  large-sized  fragments  (20  cm  and  larger)  will 
have  a  disk  shape.  The  fragment  mass  density  distribution  will  vary,  depending 
on  the  colliding  objects. 

6.2.2  Determination  of  Fragment  Spread  Velocity 

To  obtain  a  relationship  between  a  smaller  mass  of  fragment  and  its 
velocity  to  a  larger  mass  of  fragment  and  its  velocity,  it  has  been  assumed 
that  the  kinetic  energies  are  equally  imparted  such  that 

12  12 
2  '"I’^l  ■  2  '"2''2 

or 

”1  ■  Vz/ 

where  mj^  and  m2  are  arbitrary  debris  particles  with  velocities  of  v^  and  V2, 
respectively.  This  assumption  is  based  on  the  consideration  that  when  a  shock 
wave  travels,  it  applies  equal  pressure  along  the  front  at  a  given  time.  It 
causes  the  material  to  break  up  in  unequal  sizes  but  imparts  equal  kinetic 
energy.  However,  this  assumption  is  not  valid  for  smaller-sized  fragments 
(0.1  cm  and  smaller),  because  as  the  area  on  which  pressure  is  applied  by  the 
shockwave  becomes  very  small,  the  effectiveness  of  the  shockwave  is  reduced. 

For  large-sized  fragments,  this  assumption  is  in  agreement  with  the  test 
results  performed  by  PSI  and  produced  by  NASA  (Fig.  3^).  Three  curves  are 
shown  in  this  plot.  One  of  them  is  for  the  nominal  case,  and  the  other  two 
curves  are  for  extreme  cases  showing  the  upper  and  lower  limits  of  the  spread 
velocities  of  the  fragments.  In  general,  the  nominal  behavior  can  be  assumed. 


(95a) 


(95b) 
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F 

V  =  fragment  velocity 

=  relative  velocity  at  impact 


Figure  34.  Model:  On-Orbit  Debris  Total  Velocity  for 
Body-to-Body  Impact 
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It  is  also  important  to  account  for  the  total  mass  involved  before  and  after 
the  collision.  This  methodology  can  be  used  to  determine  the  number  of 
fragments,  their  masses,  and  corresponding  velocities.  Size  is  determined 
using  mass  density  criterion  for  a  fragment  of  given  mass  and  shape. 

6.3  SUMMARY 

The  methodology  presented  in  this  analysis  can  be  used  in  determining 
the  number  of  fragments,  their  masses,  and  corresponding  spread  velocities  as 
well  as  sizes.  This,  in  turn,  allows  a  better  estimate  of  the  collision 
probability  posed  by  these  debris  particles. 

6.4  PROGRAM  IMPACT 

The  algorithms  presented  in  Sections  4,  5,  and  6  have  been  coded  into  a 
user-friendly,  menu-driven  IBM  PC  software  package,  entitled  IMPACT.  This  pro 
gratr.  guides  the  user  through  the  breakup  phase  of  several  different  scenarios. 
It  allows  the  user  to  specify  such  things  as  energy  dissipated  in  a  breakup 
(percent  of  kinetic  energy  converted  to  heat  and  light),  initial  conditions  in 
several  forms  of  orbital  elements,  and  masses  of  objects  breaking  up.  Through 
out  the  duration  of  the  execution,  the  user  is  allowed  to  view  tabular  and 
graphical  outputs,  such  as  the  relationships  between  numbers  of  particles, 
sizes,  spread  velocities,  and  particle  orbital  parameters,  which  are  very 
useful  to  the  analyst.  The  resulting  output  of  this  program  is  a  properly 
formatted  input  file  for  the  program  DEBRIS  (to  be  discussed  in  Section  8). 

The  utility  of  this  software  is  that  numerous  calculations,  previously  done  on 
an  ad  hoc  basis,  are  combined  into  an  integrated,  easy  to  use,  software 
package . 
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7.  SPACE  VEHICLE  BREAKUP  DYNAMICS  IN  HYPERVELOCITY 
COLLISION— A  KINEMATIC  MODEL 

Development  of  a  hazard  model  for  on-orbit  breakups  must  necessarily  be 
approached  with  great  discipline  and  care,  since  abstractness  of  the  results 
and  difficulty  in  gathering  physical  evidence  make  it  likely  that  very  gross 
errors  in  the  model  may  pass  undetected,  even  after  the  event.  This  is 
particularly  true  in  the  case  of  near-term  hazards,  i.e.,  those  extant  during 
the  first  few  hours  and  days  after  the  breakup.  The  transient  hazards  in 
localized  regions  (near  the  "pinch  points")  typically  run  up  to  several  orders 
of  magnitude  above  the  long-term  average.  Furthermore,  their  magnitude,  loca¬ 
tion,  and  duration  are  highly  sensitive  to  the  initial  characterization  of  the 
breakup  and  the  number  and  distribution  of  initial  fragment  state  vectors 
resulting  from  the  breakup. 

Sections  5  and  6  provide,  in  principle,  a  representation  of  a  very 
simple  breakup  model  originally  developed  by  NASA  which,  until  now,  has  been 
used  extensively  in  the  orbital  debris  community.  (We  shall  therefore  refer 
to  this  subsequently  as  the  "NASA"  model,  although  this  is  an  Aerospace  inter¬ 
pretation.)  This  model  is  based  on  an  assumption  that  collision  results  are 
adequately  represented  by  a  spherically  symmetric  distribution  of  fragment 
velocities  about  the  system  center  of  mass  velocity  of  the  two  colliding 
objects.  An  inverse  relation  is  assumed  to  exist  (Fig.  3^)  between  fragment 
size  and  velocity  relative  to  the  center  of  mass.  However,  experimental  obser¬ 
vations  have  exhibited,  in  general,  a  more  complex  structure  to  the  fragment 
cloud  formation  than  this  representation  is  able  to  reflect.  Furthermore,  in 
many  cases,  large  fragments  would  have  to  survive  enormous  accelerations  to 
arrive  at  velocities  near  the  center  of  mass,  a  behavior  which  is  not  supported 
in  theory  or  observation. 

These  differences  became  important  in  the  characterization  of  hyper¬ 
velocity  collisions  of  extended  bodies,  such  as  spacecraft,  and  the  calculation 
of  the  resultant  hazards.  While  these  differences  are  largely  washed  out  in 
the  very  long-term  projections  of  fragment  populations,  overlooking  details  of 
the  initial  conditions  of  breakup  can  lead  to  significant  understatement  and/or 
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mislocation  of  peak  fragment  densities  and  resultant  hazards  in  the  near  term. 
In  this  section,  we  present  a  new  collision  model,  referred  to  as  the 
Aerospace  kinematic  model,  which  provides  an  improved  characterization  of  such 
collisions.  We  further  explore  results  of  incorporating  the  kinematic  model 
in  a  cohesive  hazard  model  for  such  events  and  contrast  these  with  results 
obtained  using  the  older  model. 

7 . 1  BACKGROUND 


Hypervelocity  collisions  of  space  vehicles  present  a  problem  whose 
study  has  been  limited  by  the  difficulty  in  accurately  modeling  such  events. 
Elaborate  hydrocode  models  are  required  for  direct  analysis  of  collisions 
between  even  the  simplest  of  objects;  the  results  generally  have  large  uncer¬ 
tainties.  When  the  colliding  objects  are  extended  and  complex  in  structure 
and  mass  properties,  such  models  tend  to  become  unmanageable.  However,  if  one 
is  interested  only  in  the  far-field  effects  of  the  collision,  some  insight  may 
be  obtained  by  examining  the  collision  dynamics  qualitatively  to  identify 
degrees  of  freedom  and  gross  constraints  on  the  mechanical  properties  of  the 
collision.  The  problem  is  thereby  simplified  to  one  of  managing  the  system's 
available  momentum,  energy,  and  mass  within  identifiable  constraints.  One  can 
then  evaluate  sensitivity  of  the  post-collision  states  to  variation  of  para¬ 
meters  associated  with  the  assumed  degrees  of  freedom.  To  accomplish  this, 
a  transfer  function  must  be  defined  relating  the  pre-  and  post-collision  states 
in  terms  of  these  parameters  as  constrained  by  the  conservation  laws.  In  this 
way,  the  range  of  possible  outcomes  of  an  event  can  be  examined  without  the 
necessity  of  modeling  each  condition  individually,  and  parametric  values  can  be 
sampled  randomly  for  statistical  analysis  of  the  system.  If  the  model  param¬ 
eters  have  been  carefully  chosen  to  have  clear  physical  significance,  then 
improved  knowledge  of  the  collision  dynamics  can  easily  be  incorporated  in  the 
form  of  revised  limits. 

Consider  the  general  properties  of  a  collision  depicted  schematically 
in  Figure  35.  The  volume  of  the  larger  object  (target)  is  divided  into  two 
regions,  shown  in  simplified  cross  section.  The  central  region  (unshaded)  is 
the  collision  volume;  i.e.,  the  volume  of  material  which  participates  directly 
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in  the  collision.  The  shaded  region  is  the  surrounding  or  noninvolved  target 
volume.  At  initial  contact,  a  set  of  hemispherical  shock  waves  is  generated 
over  the  collision  interface,  radiating  in  all  directions  from  the  locus  of 
contact  points.  As  long  as  collision  velocities  exceed  the  elastic  wave 
propagation  velocity  in  the  medium,  the  advancing  material  overtakes  the  shock 
front  ahead  of  it,  and  the  materials  of  both  objects  are  subjected  to  intense 
pressures  at  their  interface.  The  local  collision  phenomena  at  this  point  are 
exceedingly  complex  and  are  dependent  on  the  material  properties,  structure, 
and  dynamics  involved.  However,  certain  gross  properties  seem  to  be  common  to 
all  such  collisions,  strongly  influencing  the  far-field  distribution  of 
resultant  fragments.  Local  and  near-field  phenomena  are  then  important  only 
insofar  as  they  are  capable  of  shaping  the  far-field  results. 


Figure  35.  Schematic  Depiction  of  a  Collision 

The  material  at  the  interface  behaves  as  a  fluid  under  pressure;  this 
fluidity  has  the  effect  of  disassociating  the  collision  volume  and  the  sur¬ 
rounding  material  of  the  larger  object.  Nearly  all  of  the  momentum  transfer 
is  confined  to  the  collision  volume,  the  material  directly  in  the  path  of  the 
expanding  collision  front,  where  the  collision  is  totally  inelastic  and  the 
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two  objects  combine  to  behave  as  a  single  mass.  Energy  is  transferred  to  the 
noninvolved  volume  (i.e.,  the  material  not  directly  in  the  path  of  the  colli¬ 
sion)  in  a  secondary  interaction  that  involves  shock  waves,  explosive  forces, 
multiple  secondary  collisions,  and  shear  forces  in  a  transition  region  between 
the  two  zones . 

This  process  continues  until  the  collision  mass  exits  the  larger 
object,  at  which  point  the  internal  energy  of  the  collision  is  released  in  a 
radial  expansion  of  the  now  unconfined  material  in  finely  divided  fragments  in 
the  solid,  liquid,  and  gaseous  states.  The  noninvolved  volume  experiences  a 
lower  intensity,  explosive  breakup  into  fragments  of  various  sizes,  which  are 
accelerated  radially  from  its  center  of  mass  and  the  collision  axis  by  the 
fraction  of  collision  energy  transferred  to  it.  This  fraction  increases  with 
the  ratio  of  target  depth  to  projectile  dimensions  along  the  collision  axis. 

If  this  ratio  is  small,  the  target  is  fully  penetrated,  and  only  a  small  frac¬ 
tion  of  the  collision  energy  goes  into  explosive  breakup  of  the  surrounding 
target  volume. 

With  increasing  penetration  depth,  momentum  transfer  reduces  the 
relative  velocity  of  the  collision  front.  If  the  depth  ratio  is  sufficiently 
large,  the  collision  velocity  falls  below  the  elastic  wave  velocity,  and  the 
collision  is  fully  absorbed  by  the  target.  For  much  larger  targets,  ejection 
of  material  from  the  collision  takes  the  form  of  cratering,  with  the  collision 
mass  ejected  back  along  the  negative  collision  axis  and  a  net  elastic  reaction 
between  the  target  and  ejecta. 

Empirical  demonstration  of  these  properties  may  be  found  in  a  great 
variety  of  test  cases.  Figure  36  clearly  shows  (a)  the  differentiation  of  the 
collision  and  surrounding  target  volumes  in  a  small-scale,  hypervelocity 
penetration;  and  (b)  their  contrasting  breakup  characteristics  (that  of  the 
collision  volume  being  far  more  energetic).  Figure  37  shows  an  example  from 
the  opposite  end  of  the  scale,  a  head-on  collision  between  a  full-sized 
aircraft  and  missile  observed  at  White  Sands  Missile  Range  in  1978.  The 


82 


1 


"Ef'C'r'f  AfT.R  Ai-’tP 

'V(=,JC'  ^MPtTT 


Figure  36.  Flash  X-Ray  Series,  1-g  Titanium  Disc  (impacting  at 

~  4.6  km/sec  on  0.100-in.  2014-T6  aluminum  target  plate)* 


Figure  37.  Aircraft-Missile  Test,  1978  (closing  velocity  2  km/sec) 


*Photograph  from  Ref.  17 
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relative  velocity  of  the  collision  was  about  2  km/sec,  and  the  aircraf t/missile 
mass  ratio  was  about  30:1.  The  missile  was  inert  at  the  time  of  the  collision, 
since  it  had  exhausted  its  fuel  and  carried  no  warhead.  While  a  great  many 
subscale  hypervelocity  collision  tests  have  been  conducted  and  extensively 
reported  in  the  literature,  empirical  data  on  large-scale  collisions  are  very 
meager.  The  best-documented  examples  are  ground  tests  conducted  at  Arnold 
Engineering  Development  Center  by  PSI  (Ref.  13);  the  ASAT  intercept  test  (Ref. 
14);  a  HOE  flight  test  (Ref.  15);  and  the  Delta-180  flight  test  (Ref.  16). 
Because  of  instrumentation  limitations,  each  of  these  provided  only  very 
limited  data  with  respect  to  fragment  distribution. 

Qualitative  features  common  to  these  examples,  and  clearly  evident  in 
Figures  36  and  37,  include  a  distribution  of  fragment  sizes  and  velocity 
components  along  the  collision  axis  (relative  velocity  vector)  which  is  very 
nonsymmetric  and  dissimilar  to  the  distribution  of  normal  components;  the 
components  exhibit  some  degree  of  rough  cylindrical  symmetry  about  the  colli¬ 
sion  axis.  The  largest  fragments  continue  close  to  the  trajectory  of  their 
parent  object  (as  one  would  expect,  since  the  fact  that  they  remain  intact 
indicates  these  are  the  least  disturbed  by  the  collision  forces).  The  frag¬ 
ments  basically  are  organized  into  two  groups:  a  rapidly  expanding  radial 
distribution  of  large  numbers  of  finely  divided  fragments  about  the  centroid 
velocity  of  the  collision  mass,  and  a  less  energetic  expansion  of  smaller 
numbers  of  larger  fragments  about  a  centroid  close  to  the  target  velocity. 

This  suggests  a  collision  model  which  produces  a  bimodal  fragment  distribution 
incorporating  these  features  and  symmetries. 

7.2  KINEMATIC  MODEL 


In  an  inertial  reference  frame,  two  objects  of  mass  M^  and  M 

.  .  ,  .  ^ 
approach  their  joint  center  of  mass  at  velocities  and  respectively,  as 

illustrated  in  Figure  38.  It  is  convenient  to  describe  the  motion  of  the 

system  of  particles  relative  to  the  center  of  mass  (CM)  of  the  system. 
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Taking  advantage  of  S5nnmetry,  the  CM  coordinate  system  is  chosen  with  one 
axis  oriented  along  the  relative  velocity  vector,  ”  '^l  ”  ^2*  ^ 

velocity  of  Mj^  and  V  that  of  in  the  transformed  system.  Then 
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In  CM  coordinates,  the  total  linear  momentum  of  the  system  is  zero 

I  p.  =  I  m^V.  =  M^U  +  M^V  =  0  (98) 

i  i 

and  the  total  system  kinetic  energy  is 

“sys  *  I  I  "l  I  ‘  I  <”1“^  * 

Since  the  coordinate  system  is  inertial,  all  velocities  prior  to  and  following 
the  collision  are  constants.  Accelerations  during  the  collision  are  assiuned 
to  be  instantaneous. 


Immediately  after  the  collision,  the  system  is  described  as  consisting 
of  two  translating,  nonrotating,  uniformly  expanding  spheres,  each  comprising 
a  constant,  continuous,  undifferentiated  distribution  of  mass  in  velocity 
space.  These  spheres  represent  the  respective  mass  distributions  of  the 
noninvolved  target  volume  and  the  collision  volume. 

Let  the  quantities  (M|^,  U',  u,)  and  (M^,  V',  v)  identify  the  mass 
centroid  translational  velocity  and  surface  radial  expansion  rate  (spread 
velocity)  of  each  of  the  two  resultant  spheres  in  the  CM  system  (Fig.  38c).* 
Application  of  the  conservation  laws  provides  the  following  constraints. 


*In  this  section,  the  barred  notation  (u,  v,  r)  is  used  to  denote  a  specified 
value  (e.g.,  extreme  or  mean)  of  a  scalar  quantity,  rather  than  a  vector. 
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n 


Conservation  of  mass: 

M'  +  +  M2 

Conservation  of  linear  momentum; 

M^  U’  +  V’  =0 

Conservation  of  energy: 

12  2  12  12 
-  (Mj^U^  +  M2V  )  =  I  M'U'  +  I  M'V 

+  ^  J*//  pj^(u  •  u)^d^u 

Iprr  t"*  - 

+  «  iJJ  Po(v  •  v)  d  V  +  Q 
2  '^2  ^loss 


(100) 


(101) 


(102) 


In  Eq.  (102),  u  and  v  are  velocities  of  the  volume  elements  d  u  and 
3 

d  V  relative  to  the  sphere  centroids  U’  and  V,  respectively.  The  p,  (u)  and 

p2(v)  are  mass  densities  of  each  sphere  as  a  function  of  velocity  relative 

to  the  respective  sphere  centroids,  while  accounts  for  kinetic  energy 

dissipated  in  the  collision,  i.e.,  converted  to  other  forms  (heat,  light, 

etc.).  Thus,  the  first  two  terms  on  the  right  represent  the  translational 

kinetic  energy  of  sphere  centroids,  and  the  second  two  terms 

represent  the  respective  kinetic  energies  of  spread  about  the  sphere  centroids 

(KE’  ,): 

spread 


KE  =  KE’  +  KE'  .  +  Q, 

sys  trans  spread  loss 


(103) 


Definition  of  Pj^(u)  and  P2(v)  presents  some  difficulty,  since  no 
experiments  have  been  conducted  which  measure  distribution  of  material  density 
in  an  expanding  debris  cloud.  Some  experiments  (e.g..  Ref.  13)  report  observa¬ 
tions  of  material  distributed  throughout  the  expanding  cloud  volume,  but  with 
no  qviantitative  measurements  that  would  lead  to  any  conclusions  as  to  the 
specific  characteristics  of  this  distribution.  For  the  present,  therefore,  we 
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shall  assume  a  uniform  statistical  distribution  of  mass  density  over  expansion 
velocity  in  the  process  of  cloud  formation.  This  assumption  is  unlikely  to 
result  in  drastic  errors,  even  if  the  "correct"  distribution  should  prove  to 
be  different. 


(u)  =  3  M|^/4ir|ul^  =  p^ 

(104) 

(v)  =  3  M^/4ir|7|^  =  Pz 

(105) 

Conservation  of  angular  momentum  dictates  that  the  spheres  be 

nonrotating,  and  KE'  .  may  be  rewritten 

spread  ^ 


KE 


spread 


4  k 

=  2ir  I  p^u  du  +  2ir  I  p^v  dv 


3  „,-2  3  „,-2 

10  l“  10  2 


(106) 


Then  the  system  energy  is 


\  (M^U^  +  =  I  M'(U'^  +  I  i  M’(V'^  +  I  v)^  +  (107) 


7.3  TRANSFER  FUNCTIONS 

We  now  define  a  system  of  transfer  functions  expressing  the  post¬ 
collision  state  (M^,  M^,  U',V',  u,  v)  in  terms  of  the  precollision  state  (M^^, 
M^,  U,  V)  and  a  set  of  independent  parameters  governing  the  system  degrees  of 
freedom.  Assume  that  is  the  target  mass  and  is  the  projectile  mass. 
Define  a  parameter  y  as  the  fraction  of  Mj^  directly  involved  in  the  colli¬ 
sion  (0  <  Y  <  !)•  The  mass  of  the  noninvolved  target  volume  is 

m|  =  (1  -  y)M^  (108) 

Using  Eq.  (100),  the  mass  of  the  collision  volume  is  then 

M2  =  Mz  +  Y  Ml  (109) 
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Combining  Eqs .  (108)  and  (109)  with  Eq.  (101)  places  upper  limits  on  magni¬ 
tudes  of  the  resultant  centroid  translational  velocities  in  the  CM  system 


|U’I  1  |u| 


(110) 


(1  -  y)  M,  ^ 

- -  1 V 


(111) 


These  upper  limits  are  achieved  if,  in  the  collision,  no  translational  kinetic 
energy  is  converted  to  kinetic  energy  of  expansion. 


Now,  define  a  momentum  transfer  parameter  c  such  that  1  -  c  is  the 
fractional  part  of  positive  and  negative  linear  momentum  components  exchanged 
in  the  conversion  of  translational  kinetic  energy  to  kinetic  energy  of  expan¬ 
sion  (0<e  <  1  -y)*  Then 


jU'  I  =  e|S|/(l  -  y)  =  a1u| 

1^’ I  =  c|^|/[l  -  y(Mi/M2)]  =  Bjul 


(112) 

(113) 


with  A  =  e/(l  -  y)  and  B  =  e/[l  +  y(Mi/M2)]. 


The  conversion  of  translational  kinetic  energy  to  kinetic  energy  of 
expansion  must  be  apportioned  between  the  two  resultant  masses.  To  retain 
maximum  freedom  in  the  apportionment,  two  more  parameters,  a  and  (3,  are  intro¬ 
duced.  These  parameters  govern  the  kinetic  energies  of  expansion  of  the 
resultant  masses  M|  and  M^,  respectively.  Combining  Eqs.  (108)  through  (113) 
with  Eq.  (107)  then  gives  the  following  expressions  for  the  spread  velocities 
u  and  v,  respectively,  of  M^  and 

u  =  a(l  -  A)|^l  (114) 


V  =  (3(1  -  B)1V| 


(115) 
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where 


2  ^  A 
a  <  - 


M^d  -  A) 


2  3e 


+  ”2^)]  ~ 


(116) 


3^  <  - - 


Mjd  -  B) 


2  3e 


r  1  2 

-  e(M^B  +  M^A)  -  ~^- 


(117) 


In  Eqs.  (116)  and  (117),  the  equality  holds  when  mechanical  energy  is 
fully  conserved,  i.e.,  when  Qj^gg  =0  in  Eq.  (107).  In  this  case,  a  and 
3  clearly  are  not  independent;  the  value  of  one  completely  determines  the 
value  of  the  other.  This  is  true  for  any  fixed  value  of  Qj^^gg* 

Further  definition  of  the  new  parameters  a  and  (3  is  now  obtained  by 
accounting  for  Qj^^gg  and  then  ascribing  a  balance  to  a  and  |3,  which  are 
mutually  dependent. 

Considering  Eq.  (103),  a  new  parameter  n  is  defined  as  the  fraction 
of  energy  theoretically  available  to  *®’gpj.ggjj  which  is  actually  delivered 
to  ^'spread  absorbed  by  Q^^gg)*  From  Eq.  (103) 


n  =  1  - 


KE  -  ke; 
sys  trans 


KE' 

spread 


KE  -  KE’ 
sys  trans 


(118) 


Combining  this  with  Eqs.  (106),  (116),  and  (115) 


KE' 

spread 


=  yI  [M{a^  (1  -  A)^  +  M’  3^  (1  -  B)^  ] 


=  n(KE  -  KE'  ) 
sys  trans 


(119) 


Thus,  Q,  is  accounted  for  in  simultaneously  scaling  a  and  3  by  a  factor 
1/2 

n  .  Clearly,  n  always  has  a  value  between  0  and  1.  Equation  (116)  can 
be  rewritten  as  an  equality 


M^d  -  Ay 


[k^  .  -  £(MjB  .  MjA)]-  ii- 


-  „  „2 

B  ”1*^ 


(120) 
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with 


3^  < 


_5qB_ 


3eM^(l  -  B)‘ 


Ml  +  1^2  -  e(MiB  +  M2A) 


(121) 


Finally,  the  values  of  the  expansion  parameters  a  and  3  must  be 
balanced.  It  has  been  assumed  throughout  that  the  incident  bodies  are  inert, 
containing  no  significant  energy  sources  such  as  explosives  or  fuel.  (Even  if 
such  sources  are  present,  they  may  well  be  negligible  in  comparison  to  the  col¬ 
lision  energy.)  In  this  idealized  system,  then,  all  of  the  energy  of  expansion 

(KE'  ,)  for  both  the  collision  volume  (M*)  and  surrounding  volume  (M' )  must 

spread  Z  1 

originate  in  the  collision  volume.  A  proportion  of  this  energy  is  transferred 
to  the  surrounding  volume  during  the  collision  to  account  for  its  expansion, 
while  the  remainder  accounts  for  expansion  of  the  collision  volume.  (This  is 
not  intended  as  a  discussion  of  the  local  mechanics  of  interaction  between  the 
collision  volume  and  surrounding  target  volume.  Rather,  it  is  a  way  to  describe 
the  net  effects  of  the  interaction  significant  to  the  resultant  far-field  dis¬ 
tribution  of  materials.)  If  no  energy  is  transferred,  a  *  0  and  3  is  maximized; 
likewise  for  u  and  v,  respectively. 


Another  parameter,  X,  is  now  defined  as  the  fraction  of  KE'spread  trans¬ 
ferred  to  (0  <  X  <  1).  Then  the  inequality  (121)  is  replaced  by  an  equality 


Ij2  ^  -  .5n(l.-  X)B  ^  +  M  -  E  (M-B  +  M„A)] 

3eM^(l  -  zy  ^  ^  ^ 


(122) 


Note  that  q  and  X  now  have  effectively  supplanted  a  and  3  as  independent 
variable  parameters. 

Equations  (108),  (109),  (112)  through  (115),  (120),  and  (122)  form 
a  system  of  transfer  functions  which  completely  define  the  post-collision  state 
(M^,  M^,  U',  V,  u,  v)  in  terms  of  the  precollision  state  (M^,  M2,  U,  V)  and 
four  independent  parameters  (y,  c,  X,  n)*  These  parameters  form  a  simple 
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expression  of  the  system  degrees  of  freedom  in  terras  wherein  physical  signifi¬ 
cance  is  readily  grasped: 


0  <  Y  <  1 

Mass  transfer 

0  <  e  <  1-  Y 

Momentum  transfer 

0  <  X  <  1 

Energy  transfer 

0  <  q  <  1 

Energy  dissipation 

The  cloud  centroid  velocities  in  the  original  inertial  reference  frame 
are  then  given  by 


' 


U  V 


rel 


-» 

V 


cm 


V  V 


rel 


-> 

V 


cm 


u'Vrei  +  +  M^) 

'"'^rel  ^  ^  ^  ”2^ 


(123) 

(124) 


Assignment  of  values  to  the  parameters  (y*  e»  X,  n)  for  a  particular 
problem  is  dependent  upon  the  specific  geometries  and  materials  involved  and 
the  user's  assumptions  concerning  the  physics  of  the  collision.  Because  the 
physics  of  collisions  of  complex  bodies  are  poorly  understood  and  difficult  to 
predict,  the  best  approach  might  take  the  form  of  a  statistical  or  sensitivity 
analysis,  to  which  the  system  described  here  lends  itself  particularly  well. 

Application  of  this  model  is  not  restricted  to  impacts  of  a  small 
projectile  in  a  large  target;  it  is  adaptable  to  other  configurations.  For 
example,  a  glancing  or  partial  impact  of  two  bodies  of  comparable  dimensions 
may  be  handled  by  executing  twice,  allowing  each  body  in  turn  to  act  as  a 
target,  while  the  involved  portion  of  the  other  body  acts  as  a  projectile. 

7.4  CALCULATION  OF  STATE  VECTORS 


The  transition  from  a  continuous  mass  distribution  to  discrete  fragment 
states  is  accomplished  in  a  series  of  steps.  First,  the  expanding  spheres  are 
quantized  into  a  large  niunber  of  cells  of  equal  mass.  These  are  assigned 
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velocity  vectors  pseudo randomly  such  that  the  uniform,  continuous  distribution 
of  mass  over  velocity  is  replaced  by  identical,  discrete  masses  distributed 
randomly  with  uniform  probability. 


Consider  a  uniform  sphere  of  radius  R  divided  into  N  cells  of  equal 
volume.  The  volume  of  each  cell  is 

W  =  4TrR^/3N  (125) 


Assume  that  the  cells  are  subdivisions  of  M  shells  of  outer  radius 

r  (m  =  1 ,  2,  .  .  . ,  M:  r  =  0:  r..  =  R)  each  containing  n  cells.  Then 
m  o  M  °  m 

W  =  4-tr(r^  -  r^  )/3n 
m  m-1  m 


r 

m 


1/3 


The  weighted  mean  radius  of  a  shell  of  outer  radius  r  is  then 
“  m 


(126) 

(127) 


r 

m 


+  r 


m-1 


1/3 


(128) 


For  the  m-th  shell,  choose  n  sets  of  direction  cosines  (a.  ,  b.  ,  c.  )  as 

m  111 

m  m  m 

uniformly  distributed  random  numbers  between  -1  and  1,  and  their  opposites 

(-a.  ,  -b.  ,  -c .  ),  for  i  =1,  2,  ...»  n  /2.  The  uniform  sphere  is  then 
1  1  1  m  m 

m  m  m 

approximated  by  the  N  cells  of  volume  W  having  position  vectors 


i  c^(a.  i,  b.  j,  c.  k). 
m  m  m 


A  cell  mass  m  is  selected;  the  sphere  representing  the  noninvolved 
volume  is  then  approximated  by  setting  N  =  M^'/m,  R  =  u.  For  the  collision 
volume,  set  N  =  M^'/m,  R  =  v.  This  method  absolutely  conserves  momentum  and 
approximately  maintains  the  mass  distribution.  Mass  distribution  and  energy 
content  cannot  be  maintained  simultaneously  in  converting  from  a  continuous  to 
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a  discrete  and  finite  distribution  of  masses;  therefore,  this  process  necessi¬ 
tates  a  small  loss  of  system  energy.  This  generally  may  be  accounted  for  by 
incorporation  in 

Next,  a  list  of  discrete  fragment  masses  is  generated  for  each  volume 
(i.e.,  the  collision  volume  and  the  noninvolved  volume).  The  approach  employed 
here  is  based  on  relationships  documented  in  Reference  13.  However,  in  incor¬ 
porating  this  approach  in  the  kinematic  model,  we  have  further  developed  and 
extended  these  relationships  to  form  a  self-consistent  method  of  calculating  a 
fragment  list  for  any  hypervelocity  collision.  In  the  process,  potential 
inconsistencies  observed  in  previous  applications  of  these  relationships  are 
eliminated. 

Reference  18  identifies  differing  relationships  for  fragment  mass 
distributions  resulting  from  explosive  breakup  and  collision  fragmentation  as 
follows . 


Explosion  Fragments; 


N  (m  >  ra^)  =  KjM  exp 


(129) 


where  N  is  the  cumulative  number  of  fragments  having  mass  larger  than  m^,  and 
M  is  the  total  mass  of  the  fragmented  object.  The  constants  and  are 
defined  differently  for  the  largest  fragments  and  those  in  a  smaller  size  range, 
based  on  empirical  data. 


>  1936 

gm: 

=  1.71 

X 

o 

1 

=  0.02056 

<  1936 

gm: 

=  8.69 

X  10 

=  0.05756 

Collision  Fragments; 

(m  >  m  )  =  K(M  /m  )P  (130) 

o  e  o 

where  K  and  p  are  constants  defined  empirically  as  K  =*  0.4478,  p  =  0.7496. 
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Here,  is  the  "ejected  mass,"  i.e.,  the  mass  which  would  be  excavated  from 
the  target  in  a  cratering  impact  if  the  target  were  very  much  larger  than  the 
projectile  striking  it,  defined  as 


M  =  V^m  (131) 

e  p 

where  m^  is  the  projectile  mass,  and  V  is  the  relative  velocity  magnitude  at 
impact  in  Km/s- 

As  stated  earlier,  some  inconsistencies  have  been  found  in  past  applica¬ 
tions  of  Eqs.  (127)  through  (129)  to  hypervelocity  collision  problems.  Ref¬ 
erence  19  notes  that  the  larger  fragments  of  a  collision  event  are  distributed 
per  Eq.  (129),  while  the  smaller  fragments  are  distributed  per  Eq.  (130),  so 
that  the  total  distribution  is  a  combination  of  the  two.  However,  no  rationale 

is  given  as  to  how  the  two  distributions  should  be  combined.  Also  problematic 

2 

is  the  fact  that  unless  tho  target  has  at  least  V  times  the  mass  of  the 
projectile,  Eqs.  (130)  and  (131)  produce  more  than  the  available  mass  in 
fragments.  Reference  18  circumvents  this  by  stating  that  the  target  mass  must 
be  greater  than  10  x  M^,  but  Reference  19  ignores  this  condition,  applying  the 
distribution  to  a  breakup  where  the  target  and  projectile  are  of  comparable 
mass.  Our  integration  of  these  observations  with  the  kinematic  breakup  model 
incorporates  a  proposed  resolution  of  these  inconsistencies. 


In  the  kinematic  model,  the  low-intensity  explosive  breakup  Eq.  (129) 

is  applied  to  the  noninvolved  target  mass  To  generate  a  fragment  list,  a 

minimal  fragment  size  m^  is  identified,  and  the  total  number  of  fragments  in 

a  range  between  m.  and  m.  ,  is  identified  as  follows: 

1  1+1 


n(m.  <  m 
1 


<  m.  , )  =  K,M 
—  1+1  1 


,  „  ,1/2  1/2.  X  V  f  1 ^l/2  1/2.. 

exp(-  K  1  m  )  -  exp[-  K_(i  +1)  m  )] 

+  O  L  O 


=  n  (mp  (132) 

where  m^^  =  i  m^;  mj^  =  (2i  +  l)m^/2.  A  list  is  generated  with  M  =  using 
each  set  of  constants  (K^^,  K^).  Fragments  larger  than  1936  gm  are  kept  from 
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the  (m  >  1936  gm)  list,  and  fragments  from  the  (m  <  1936  gm)  list  are  kept 
from  1936  gm  downward  until  ^  (m  <  1936  gm)  =  J  (m  >  1936  gm).  The  two 

lists  are  then  combined. 


The  collision  fragment  distribution  Eq.  (130)  is  similarly  applied  to 
the  collision  mass  M^.  A  list  is  generated  by 

n(m:)  =  K(vJ^.M,/m^)P  [i'P  -(i  +  1)'^]  (133) 

1  rei  z  o 

for  very  small  m^.  Then  fragment  masses  are  accumulated  from  i  =  0  up  to 
i  =  I  such  that 


I 

1  m:  =  Ml  (134) 

i=0 


and  the  remainder  are  discarded. 
2 


when  M  = 
e 


'^Rel  ”2  "  ”l  "  ”2- 


This  eliminates  the  problem  posed  by  Eq.  (131) 


Note  that  Eq.  (131)  is  dimensionally  incorrect.  It  is  an  empirically 

derived  relationship  which  introduces  the  effect  of  collision  energy  variation 

2 

into  the  resultant  fragment  mass  distribution.  Without  the  V  factor, 

Eq.  (130)  would  produce  the  same  fragment  distribution  for  any  projectile/ 
target  combination,  regardless  of  the  kinetic  energy  of  impact.  The  interpre¬ 
tation  employed  in  the  kinematic  model  preserves  the  effect  of  increased 
collision  energy  resulting  in  a  breakup  into  larger  numbers  of  smaller-sized 
fragments . 


The  basic  cells  must  be  no  larger  than  the  smallest  fragment  of  interest 
from  the  fragment  list.  Combinations  of  cells  are  now  accumulated  to  produce 
fragments  approximating  the  fragment  masses  on  the  fragment  list.  This  is 
accomplished  by  a  random  selection  process  which  preserves  the  expected 
inverse  relationships  between  fragment  mass  and  velocity  relative  to  the  inci¬ 
dent  velocity  of  the  body  from  which  the  fragment  originates.  Each  fragment 
on  the  fragment  list  thereby  is  assigned  a  velocity  which  is  the  vector  sum  of 
the  velocities  of  its  component  cells. 
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7.5 


FRAGMENT  DISTRIBUTION 


As  fragment  size  decreases,  the  numbers  increase  geometrically,  so  that 
tracking  the  individual  state  vectors  of  fragments  becomes  very  impractical 
when  one  is  interested  in  those  of  very  small  size.  Furthermore,  in  order  to 
utilize  the  results  for  statistical  analysis,  it  is  necessary  to  identify  a 
statistical  distribution  n(v,d),  i.e.,  fragment  number  density  as  a  function  of 
fragment  velocity  v  and  size  d.  Although  the  overall  fragment  distribution  is 
not  (in  general)  spherically  symmetric,  often  it  can  be  broken  down  into  com¬ 
ponents  which  can  be  treated  as  spherically  symmetric  about  their  own  centroids. 
Then  the  vector  velocity  v  c<in  be  replaced  by  a  velocity  magnitude  v  =  |  v  j 
relative  to  the  appropriate  centroid,  so  that  the  distribution  is  simplified  to 
a  function  n(v,d)  of  two  rather  than  four  variables.  The  overall  fragment 
distribution  can  then  be  constructed  by  the  principle  of  superposition  of  the 
component  distributions. 

The  output  of  the  kinematic  model  is  a  distribution  of  number  versus  size 
n^(d)  of  all  fragments  above  the  minimal  size  of  interest,  individual  frag¬ 
ment  velocities  down  to  the  smallest  size  which  is  computationally  manageable, 

^  mm 

and  the  centroid  velocities  (O’,  V')  and  surface  radial  expansion  rates  (u,  v) 
for  each  of  two  spherical  fragment  distributions.  The  larger  fragment  veloc¬ 
ities  can  be  analyzed  to  obtain  their  statistical  distribution  n^(v,  d);  but 
this  distribution  must  be  inferred  for  the  small  fragments,  for  which  only  the 
upperbound  velocities  (u,  v)  are  given.  Once  a  functional  form  n^(v,d^)  is 
identified  for  each  size  range  d^^  +  Ad,  the  fragments  of  each  size  category 
within  a  velocity  range  (v^,  v^)  can  be  enumerated: 

nd(di  +  Ad)  =  J  nv(v,di)dv  (135) 

VI 

Figure  39  illustrates  an  example  of  a  fragment  distribution  generated 
by  the  kinematic  model.  Velocities  are  spherically  symmetric  about  the  cen¬ 
troid  of  expansion.  A  simple  three-dimensional  approximation  of  the  distribu¬ 
tion  of  fragments  larger  than  1  mm  is  obtained  by  superposition  of  three 
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V  m/sec 


isotropically  distributed  spheres  of  particles  in  velocity  space  having  the 
following  characteristics: 


Radius 
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Figure  39.  Fragment  Distribution 
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7.6 


EFFECT  OF  SECONDARY  COLLISIONS 


So  far,  it  has  been  assumed  that  the  collision  and  breakup  occur  instan¬ 
taneously,  and  that  the  collision  fragments  emanate  radially  from  a  point 
source.  In  reality,  the  colliding  bodies  have  finite  dimensions,  and  the 
collision  occurs  over  a  finite  time  interval.  While  these  properties  can  be 
ignored  in  describing  the  general  motion  of  the  resultant  system  of  fragments, 
they  do  significantly  influence  the  detailed  distribution  of  fragments  within 
the  system.  This  is  due  to  secondary  interactions  among  the  different  com¬ 
ponents  of  the  system  very  early  (on  the  order  of  milliseconds)  after  the 
collision,  when  the  dimensions  of  the  colliding  bodies  are  still  significant 
compared  to  the  distances  traveled  by  the  fragments. 

Of  particular  importance  is  the  physical  observation  that  the  larger 
fragments  are  generated  external  to  the  volume  of  the  direct  colliding 
masses.  If  the  spread  velocity  of  the  fragments  emanating  from  the  hyper¬ 
velocity  collision  /olume  exceed  the  relative  velocity  of  the  centroids  of  the 
collision  and  explosive  breakup  volumes,  some  fraction  of  the  hypervelocity 
collision  fragments  will  overtake  and  undergo  secondary  collisions  with  the 
explosive  fragments.  Assuming  a  fair  degree  of  inelasticity  in  these  secondary 
interactions,  one  can  expect  the  collision  fragments  to  give  up  much  of  their 
kinetic  energy  (thereby  lowering  their  velocities)  relative  to  the  generally 
larger  and  much  more  massive  explosive  fragments.  This  interaction  also  will 
impart  momentum  to  the  explosive  fragments  in  directions  radial  to  the 
collision  centroid,  but  the  large  total  mass  ratio  of  the  explosive  fragments 
to  the  involved  collision  fragments  results  only  in  a  small  net  velocity 
change  to  the  larger  fragments. 

Physical  evidence  of  this  kind  of  interaction  is  found  in  an  unreported 
observation  from  the  PSI  ground  tests.  In  tests  that  had  full  penetration  of 
the  target,  resulting  in  a  collision  centroid  in  positive  relative  motion  along 
the  collision  axis  from  the  target  body,  it  was  observed  that  the  explosive 
breakup  fragments  from  the  noninvolved  portion  of  the  target  body  received  a 
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net  momentum  component  along  the  negative  collision  axis  (in  the  direction 
from  which  the  projectile  had  come),  and  large  numbers  of  collision  fragments 
were  found  in  and  around  the  larger  fragments.  While  no  attempt  was  made  to 
quantify  these  observations,  this  negative  component  in  the  motion  of  the 
large  target  fragments  is  difficult  to  explain  except  by  secondary  momentum 
transfer  from  the  spreading  collision  fragments. 

As  a  result  of  such  a  secondary  interaction,  it  may  be  expected  that 
the  distribution  of  the  involved  fraction  of  the  collision  fragments  will  be 
concentrated  in  the  region  of  the  explosive  fragments,  increasing  the  overall 
fragment  density  in  that  region.  Calculation  of  this  resultant  distribution 
can  be  separated  into  two  component  problems:  calculation  of  the  fraction  of 
collision  fragments  involved,  and  calculation  of  the  resultant  changes  in 
individual  fragment  momenta.  In  both  cases,  exact  solution  would  represent  a 
very  complex  problem  in  statistical  mechanics.  Since  precision  and  accuracy 
of  results  is  no  better  for  such  a  problem  than  that  of  the  input  information 
(fragment  numbers,  sizes,  shapes,  velocities,  etc.),  which  in  this  case  can 
only  be  grossly  estimated,  no  such  analysis  is  warranted  or  attempted  here.  A 
first  approximation  is  arrived  at  by  much  simpler  methods. 

Suppose  the  collision  phase  of  a  breakup  is  completed  at  time  t  =  0 
such  that  the  collision  fragments  are  uniformly  distributed  within  a  spherical 
volume  whose  surface  is  expanding  uniformly  at  radial  velocity  of  constant 
magnitude  R.  At  t  =  0  the  surface  radius  is  R^.  Assume  also  that  a  body 
having  apparent  cross-sectional  area  A  (as  seen  from  the  sphere  centroid)  is 
initially  at  distance  R^  from  the  sphere  center  and  moving  radially  at  constant 
speed  V  <  R  relative  to  the  sphere  centroid.  Then,  at  any  time  after  t  =  0, 
the  flux  of  particles  through  area  S  (i.e.,  colliding  with  the  external  body 
per  unit  time)  is 

(t)  =  p(t)  S  AV(t)  (136) 

where  p(t)  is  the  number  density  collision  fragments  of  per  unit  volume  of 
the  sphere,  and  AV(t)  is  the  velocity  differential  between  the  external  body 
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and  fragments  at  a  radius  +  Vt  from  the  sphere  centroid.  If  N  is  the 
total  number  of  collision  fragments 


p(t ) 


V(t) 


3N _ 

Air(R  +  Rt)^ 
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’r  +  Vt 

R  °  -  -- 

R  +  Rt 
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(137) 
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The  fraction  of  collision  fragments  involved  in  secondary  collisions  with  the 
external  body  after  t  =  0  is  then  given  by 
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=  ~~  (  1  -  M  (V  <  R)  (139) 

WR^  \  R/ 

This  assumes  that  S  is  constant.  Further  levels  of  detail  (e.g.,  time  varia¬ 
tion  of  S,  oblique  orientation  of  surfaces,  perpendicular  component  of  initial 
displacement)  are  second-order  considerations  which  do  not  greatly  alter  the 
result,  except  in  the  case  where  V  is  close  to  (or  exceeds)  R;  this  would  lead 
to  a  unique  mathematical  formulation  for  each  scenario  and  geometric  config¬ 
uration  studied. 


The  form  of  the  scattering  function  likewise  is  unknown,  being  dependent 
on  a  great  number  of  factors  such  as  size,  shape,  composition,  and  orientation 
of  the  scattering  (large)  fragments.  Elasticity  of  the  scattering  should  be  a 
function  of  the  scattering  angle;  scattering  near  180®  from  the  direction  of 
incidence  is  highly  inelastic,  while  scattering  at  small  angles  may  be  highly 
elastic.  From  a  classical  scattering  theory,  if  the  scattering  is  treated  as 
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isotropic  on  the  average,  then  a  crude  estimate  of  the  resultant  velocity 
distribution  can  be  obtained.  If  it  is  assumed  that  scattering  at  angles 
between  100®  and  180°  is  essentially  inelastic,  the  fraction  of  incident 
particles  scattered  in  this  range  is 


1 

2 


sin0d6~O.4 


(140) 


Thus,  40^  of  the  scattered  particles  will  give  up  all  momentum  relative 
to  the  scattering  fragments,  resulting  in  their  capture  within  the  velocity 
region  of  the  large  fragments.  The  remainder  are  assumed  to  be  distributed  at 
relative  velocities  out  to  R  -  V  from  the  centroid  of  the  scatterers. 

As  an  example.  Figure  40  illustrates  the  incident  and  resultant  frag¬ 
ment  velocities  (as  seen  from  above)  from  a  hypothetical  glancing  collision 
between  two  orbiting  bodies  labeled  A  and  B,  whose  orbits  are  mutually  inclined 
at  20°.  The  bodies  are  assumed  to  be  generally  cylindrical,  of  similar  dimen¬ 
sions,  but  with  a  mass  ratio  =  5/3.  The  collision  is  assumed  to  occur 

with  the  body  axes  oriented  parallel  to  the  relative  velocity  vector  and 
involve  =  Yg  =  1/3  the  mass  of  each  body.  Figure  41  details  the  relative 
positions  of  the  collision  fragments  and  those  of  the  lower  intensity  explo¬ 
sive  breakup  of  the  noninvolved  portion  of  object  A  during  the  first  second 
after  the  collision.  If,  in  Eq.  (139) 

R  =  1  m 

o 

R  =  1 . 7  x  10^  m/s 

V  =  10^  m/s 

S  =  3  m^ 

then  the  fraction  of  collision  fragments  intercepted  by  A  is  q.  »  10%.  From 

A 

Eq.  (140),  about  4%  of  the  total  might  be  distributed  within  the  velocity 
regime  of  the  large  fragments.  The  remaining  6X  is  estimated  to  lie  within 
approximately  500  km/s  of  the  centroid  of  the  large  fragments. 
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PARALLEL  TO  V^g,  (km/s) 


RELATIVE  POSITION  (m) 


Figure  41.  Early  Evolution  of  Debris  Distribut 
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TIME 


=  0.01  sec 


TIME  =  1  sec 


Figure  Al.  Early  Evolution  of  Debris  Distribution  (Cont 
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Clearly,  in  this  scenario.  Object  B  represents  a  case  where  V  »  R  and 
higher-order  effects  would  dominate, 

7.7  COMPARISON  OF  RESULTS  WITH  OTHER  MODELS 


The  kinematic  model  has  been  employed  to  provide  the  hazard  analysis 
predictions  on  two  experiments  to  date  involving  actual  spacecraft  collisions: 
an  ASAT  test  conducted  in  September  1985  and  the  Delta  180  test  in  September 
1986.  In  both  cases,  analysis  of  the  results  showed  very  good  agreement  with 
the  kinematic  model  predictions  (Refs.  14,  16,  20,  and  21). 

In  contrasting  these  results  with  predictions  based  on  the  "NASA"  model, 

some  important  differences  are  observed.  In  Section  2,  Tables  1  and  2  are  two 

predictions  of  the  results  of  the  same  collision.  Table  1  is  based  on  the 

older  "NASA"  model,  while  Table  2  was  generated  by  the  kinematic  model  and 

reflects  Figure  39,  neglecting  fragments  which  will  not  survive  one  orbit. 

Figure  42  shows  the  probability  of  collision  with  a  fragment  larger  than  1  mm 

2 

in  diameter  for  a  spacecraft  of  20  m  cross  section  passing  through  the  cen¬ 
troid  of  the  debris  cloud  represented  by  Table  2,  at  90*  mutual  inclination 
(total  angle  between  the  orbit  planes)  with  the  centroid  orbit.  These  results 
were  obtained  using  the  propagated  cloud  hazard  model  incorporated  in  the 
kinematic  model  companion  tool.  Program  DEBRIS  (described  in  the  following 
section).  Figure  42  shows  the  hazard  for  a  single  passage  through  the  debrio 
cloud  center  at  any  time  during  the  first  24  hr  after  the  collision  event. 

The  most  significant  features  of  Figure  42  are  the  very  large  spikes  in 
collision  probability  encountered  at  even  revolutions,  when  the  cloud  is 
passing  through  the  primary  "pinch  point,"  along  with  lesser  rises  observed  at 
the  half-revolution  points.  It  is  clear  that  time  averaging  of  collision 
probabilities  should  be  avoided  when  one  is  concerned  with  the  actual  hazard 
to  be  experienced  by  a  spacecraft  at  a  specific  location  and  time.  While  the 
"average"  hazard  quickly  falls  below  the  level  of  hazard  posed  by  the  equi¬ 
valent  natural  background  fltjc  of  meteroids  (dashed  line),  the  peaks  remain 
several  orders  of  magnitude  above  this. 
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ORBITAL  REVOLUTIONS  OF  CLOUD  CENTROID 


Figure  A2,  Collision  Expectation  versus  Cloud  Position 

Figure  A3  compares  the  propagated  results  based  on  initial  conditions 
generated  by  the  kinematic  model  with  those  of  the  "NASA"  model.  The  upper 
curve  (kinematic  model)  of  Figure  43  is  identical  to  Figure  42,  with  the 
horizontal  scale  expanded;  the  lower  curve  is  based  on  Table  1,  similarly 
propagated.  Two  significant  differences  are  observed: 


The  overall  lower  level  of  hazard  predicted  by  the  "NASA"  model 
during  this  interval,  despite  predicting  a  greater  number  of 
orbiting  fragments,  is  due  to  the  fact  that  all  but  a  small 
fraction  of  these  are  assigned  very  high  spread  velocities, 
resulting  in  their  more  rapid  dispersal. 

Displacement  of  the  "NASA"  model  peaks,  relative  to  those  of  the 
kinematic  model,  is  due  to  the  "NASA"  model  forcing  the  cloud  cen¬ 
troid  and  larger  fragments  to  the  system  center  of  mass  of  the  two 
colliding  bodies,  while  the  kinematic  model  possesses  additional 
degrees  of  freedom  to  permit  a  bimodal  distribution  of  multiple 
clouds  with  separate  centroids.  (In  the  particular  case  shown, 
the  second  cloud  centroid  is  suborbital,  and  its  associated 
fragments  do  not  contribute  significantly  to  the  on-orbit  hazards 
beyond  the  initial  expansion. ) 
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ORBITAL  REVOLUTIONS  OF  CLOUD  CENTROID 


Figure  43.  Comparison  of  Kinematic  Model  with  "NASA"  Model 
7.8  CONCLUSIONS 


As  stated  at  the  beginning  of  this  section,  the  calculation  of  near-term 
on-orbit  debris  hazards  resulting  from  a  collision  of  two  spacecraft  is  highly 
sensitive  to  initial  characterization  of  the  breakup.  Improvements  in  subse¬ 
quent  propagation  of  the  debris  clouds  cannot  improve  the  results  if  the 
initial  conditions  are  wrong.  The  results  shown  in  Figure  43  illustrate  the 
differences  produced  by  differing  characterizations  of  the  same  event,  both  in 
the  overall  resultant  hazard  level  and  in  the  location  of  hazard  peaks. 

Identification  of  the  "correct"  overall  or  "average"  level  of  hazard  is 
highly  debatable  and  will  remain  so  until  much  more  definitive  experimental 
data  on  large-scale  hypervelocity  collisions  become  available.  More  work  needs 
to  be  done  in  this  area;  however,  high  cost  has  been  a  factor  in  limiting 
experimentation  of  this  nature. 
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Section  7.7  shows  plainly,  however,  that  accuracy  in  locating  the  narrow 
regions  of  peak  hazard  is  of  more  importance  than  the  particular  numerical 
value  obtained.  It  is  passage  through  these  regions  that  will  pose  a  signifi¬ 
cant  hazard  to  resident  spacecraft,  regardless  of  the  precise  numerical  level 
of  that  hazard.  Protection  of  spacecraft  from  near-term  hazards  depends  on 
the  ability  to  predict  and  avoid  coincidences  with  these  regions. 

Experimental  evidence  to  date  indicates  that  the  kinematic  model  pro¬ 
vides  credible  representations  of  the  distributional  structure  (and  therefore 
peak  hazard  regions)  of  debris  clouds  resulting  from  a  spacecraft  collision. 
This  should  be  sufficient  to  justify  its  further  application  to  such  problems. 
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8.  DESCRIPTION  OF  PROGRAM  DEBRIS 


8.1  INTRODUCTION 


Program  DEBRIS  is  described  in  this  section.  The  program  determines  the 
intervals  during  which  a  given  payload  satellite  travels  through  an  expanding 
debris  cloud  and  calculates  the  probability  of  collision  associated  with  each 
transit.  This  section  provides  a  general  overview  of  program  execution.  Also 
included  are  detailed  analyses  of  the  algorithms  used  in  DEBRIS.  The  necessary 
inputs  to  use  the  program  are  discussed,  along  with  the  types  of  information 
generated  as  output. 

Program  DEBRIS  was  developed  to  calculate  the  short-  and  long-term 
probability  of  collision  of  a  given  payload  satellite  whose  orbit  is  in  the 
vicinity  of  a  newly  formed  debris  cloud.  The  program  determines  the  intervals 
during  which  the  satellite  travels  through  the  expanding  debris  cloud  and 
calculates  the  probability  of  collision  associated  with  each  transit.  These 
probabilities  are  then  combined  to  obtain  the  cumulative  probability  of 
collision,  which  is  compared  against  the  meteorite  background  to  determine  if 
a  significant  hazard  exists. 

An  unlimited  number  of  payload  satellites  may  be  examined  by  DEBRIS 
during  program  execution.  A  single  maneuvering  vehicle,  such  as  the  space 
shuttle,  may  also  be  simulated. 

The  expanding  debris  cloud  is  propagated  using  the  model  developed  in 
previous  sections.  A  number  of  sub-clouds,  each  composed  of  particles  of 
uniform  size,  distribution,  and  separation  velocity,  have  been  used  to 
represent  a  large  debris  cloud  containing  a  diversity  of  particles. 

8.2  PROPAGATION  OF  THE  DEBRIS  CLOUD 

At  time  tj,  a  target  spacecraft  breaks  up,  causing  the  formation  of  a 
debris  cloud.  At  any  subsequent  time  t,  the  debris  cloud  is  modeled  as  a 
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torus  extending  along  the  trajectory  of  the  target  orbit  with  an  elliptical 
cross  section  of  varying  dimension  (Fig.  44).  The  center  of  the  debris  cloud 
travels  along  an  orbit  with  the  same  elements  as  the  former  orbit  of  the 
spacecraft.  The  in-plane  arc  displacement  between  the  leading  and  trailing 
edges  of  the  cloud  increases  linearly  with  time  until  the  torus  closes.  The 
cross-sectional  dimensions  of  the  cloud  are  functions  of  the  total  angular 
displacement  between  the  point  of  interest  and  the  point  where  the  target 
spacecraft  disintegrated. 

PAYLOAD  SATELLITE  TRAJECTORY 


Figure  44.  Debris  Cloud  Propagation  Model 

The  volume  of  the  debris  cloud  may  be  calculated  using  either  of  two 
methods.  For  purposes  of  calculating  the  probability  of  collision  associated 
with  the  passage  of  a  payload  satellite  through  the  cloud»  the  smaller  of  the 
two  volumes  will  be  used.  This  ensures  that  the  collisional  hazard  due  to  the 
cloud  is  not  underestimated. 

It  is  possible  to  represent  a  debris  cloud  as  the  aggregate  of  several 
sub-clouds  (Fig.  45).  Each  sub-cloud  is  composed  of  a  different  number  of 
debris  particles,  of  uniform  size  and  propagation  rate.  The  superposition  of 


Figure  45.  Combination  of  Several  Debris  Clouds  to  Represent  an 
Aggregate  Cloud  of  Varying  Growth  Rates 

several  sub-clouds,  each  with  a  different  shape  and  density,  is  equivalent  to 
a  larger  cloud  of  nonuniform  density. 

8.3  EQUATIONS  FOR  DEBRIS  CLOUD  PROPAGATION  AND  VOLUME  CALCULATION 

For  any  time  t  measured  from  the  time  of  cloud  formation  tj,  the 
dimensions  and  volume  of  the  debris  cloud  may  be  calculated  using  the 
following  equations. 
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It  is  assumed  that  the  debris  cloud  grows  at  a  uniform  rate  along  the 
trajectory  of  the  target  orbit.  For  times  prior  to  torus  closure,  the 
in-plane  arc  distance  between  the  leading  and  trailing  edges  of  the  debris 
cloud  is  eq\ial  to 

L(t)  =  (L^  +  L^)t  (141) 

where 

=s  user-specified,  leading-edge  growth  rate,  Icm/min 
=  user-specified,  trailing-edge  growth  rate,  km/min 


The  circumference  C  of  an  ellipse  of  semi-major  axis  a  and  eccentricity 
e  is  given  by 


C  =  2ir  V (1/2)  a^(2  -  e^) 


The  time  at  which  the  torus  closes  (t^j^)  may  then  be  calculated  as 


(142) 


CL 


*  '■2> 


(143) 


where  C  is  the  circumference  of  the  target  orbit. 


For  all  times  t  <  t^j^,  the  length  of  the  torus  L  is  equal  to 


L 


(144) 


At  any  given  point  within  the  leading  and  trailing  edges  of  the  cloud, 
the  boundaries  of  the  elliptical  cross  section  are  defined  by  Reference  22  as 


a(t,0)  = 


(145) 


and 


b(t,e)  =  a33{t,e)  (^) 


(U6) 


where 


a(t.e) 

b(t,9) 

0 

Av 

03 

^21’^22’®33 


semi-major  axis  of  elliptical  cross  section,  km 

semi-minor  axis  of  elliptical  cross  section,  km 

total  angular  displacement  from  point  where  target  was 
broken  up  to  point  of  interest  [see  Eq.  (170)] 

expansion  velocity  of  debris  particles,  km/sec 

orbital  mean  motion  of  the  reference  orbit,  rad/sec 

functions  previously  introduced 


In  order  to  calculate  the  probability  of  collision  associated  with  the 
passage  of  a  payload  satellite  through  the  debris  cloud,  it  is  necessary  to 
determine  the  volume  of  the  cloud  at  different  times.  Assume  that  a  satellite 
enters  the  cloud  at  time  t^  and  angular  displacement  0^,  exiting  the  cloud  at 
*’N’  <  t  <  tj^  and  0^  <  0  <  tj^,  the  cloud  volume  is 


V0L^(t,e)  =  ^ 


|a3j(t,0)  a,,(t,e)|  +  [a,,(t,e)] 


22 


21 


* 


W) 


where 

aii(t,0)  =  -3(0  t  +  4  sin  0  (148) 
a2i(t,0)  =  2[1  -  cos  0  +  g3(t,0)  (1  +  cos  0)]  (149) 
a22(t,0)  =  sin  0  +  g2(t,0)  [SGN(sin  0)  -  sin  0]  (150) 
a33(t,0)  =  g3(t)  +  jsin  0|  (151) 
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and 


g^(t)  =  MINj[C^t  +  C^Cut  -  sin  8)],  ij 
g^Ct)  =  MIN  I  |c^t  +  |c^(«t  -  sin  8),  l| 
g^Ct)  =  MINjc^t,  (a  0)  sin  i  )(Av  ) 

where 

i  =  inclination  of  reference  orbit,  rad 

=  constants  discussed  in  Section  h 
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Figure  46.  Geometry  of  a  Typical  Pass  Through  the  Debris  Cloud 
The  function  SGN(X)  is  defined  as 

SGN(X)  =  1^,  if  X  0 
SGN(X)  =0,  if  X  =  0 


(152) 

(153) 

(154) 
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The  functions  gp  g^.  and  g^  model  the  effects  of  and  drag  on  cloud 
propagation.  As  previously  shown  in  Section  U,  the  constant  Cj^  incorporates 
the  effect  on  the  line  of  apsides,  shows  the  effect  on  the  line  of 
nodes,  and  represents  the  effect  of  atmospheric  drag.  A  second  method  for 
calculating  the  volume  of  the  debris  cloud  for  (t,0)  is 

VOL^Ct,©)  =  ir  a(t,0)  b(t,0)  L(t)  (155) 

where  a,  b,  and  L  are  calculated  using  Eqs.  (lAl)  and  (142). 

If  the  difference  between  the  total  angular  displacement  at  entry 
(0  )  and  exit  (0.,)  is  less  than  10®,  then  the  two  volume  functions  are 

O  N 

evaluated  at  entrance  and  exit  and  averaged 

=  i[V0L,  (t  ,  Q)  +  VOL,  (t„,  0,^)]  (156a) 

12  loo  INN 

^2  “  ®o^  V0L2(tj^,  0J,)]  (156b) 

Otherwise,  the  two  volumes  are  calculated  from 

_  1  N 

VOL^  =  S  Z  V0L^(t^,  0.)  (157a) 

i=o 

_  ,  N 

VOL^  =  J  Z  V0L2(t.,  0J  (157b) 

i=o 

where 

N  .  L  .  INT((le^  -  9^1  >  21.  INTt(t^  -  t^)(T)  ‘'l)  /  ^  ) 

h  ■  ‘o  ‘  <'n  -  'o>  (s)  (i  =  0.  1.  .  .  N) 

©i  =  determined  through  propagation  of  payload  orbit 
T  =  period  of  reference  orbit 
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The  smaller  of  VOL^^  and  VOL2  is  used  for  the  computation  of  the 

probability  of  collision  associated  with  (t  ,  t  ).  However,  if  (t  +  t  )/2 

_  o  N  ON 

is  greater  than  t  ,  then  VOL„  is  always  used  as  the  volume  of  the  debris 

^  InI  fc 

cloud . 

8.4  DETERMINATION  OF  TRANSIT  TIMES  OF  A  SATELLITE 
THROUGH  THE  DEBRIS  CLOUD 

The  time  intervals  during  which  a  satellite  is  within  the  boundaries  of 
the  debris  cloud  may  be  determined  from  the  simultaneous  solution  of  two  in¬ 
equalities  for  all  possible  values  of  t 


(t)  <  0 

(158) 

(t)  <  0 

(159) 

The  function  of  q^^Ct)  represents  the  distance  from  the  payload  satel¬ 
lite  to  the  leading  or  trailing  edge  of  the  cloud,  measured  in  the  orbital 
plane  of  the  debris  cloud  (see  Fig.  47).  The  function  q2(t)  represents  the 
out-of-plane  distance  from  the  payload  satellite  to  the  cross-sectional  bound¬ 
ary  of  the  debris  cloud  (see  Fig.  48).  It  is  necessary  that  both  of  these 
functions  be  nonpositive  for  the  satellite  to  be  inside  the  debris  cloud  at  a 
given  time. 


DEBRIS  CLOUD 


Figure  47.  Illustration  of  qj^(t) 
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PAYLOAD  SATELLITE 


Figure  48.  Illustration  of  q2(t) 


The  global  set  of  time  intervals  corresponding  to  multiple  passes  of  a 
satellite  through  the  debris  cloud  may  be  found  by  solving  Eqs.  (158)  and 
(159)  separately  and  taking  the  intersection  of  the  two  sets  of  intervals  as 
the  result.  This  is  accomplished  in  Program  DEBRIS  using  a  Newton-Raphson 
iterative  method.  The  derivatives  of  q^^  and  q^  are  numerically  evaluated, 
with  the  convergence  tolerance  set  to  10~^  km. 

For  any  given  time  t,  the  position  (^)  and  velocity  (R)  of  the  center 
of  the  debris  cloud  can  be  determined  by  propagating  the  position  and  velocity 
of  the  target  spacecraft  from  time  t^  to  t.  A  unit-vector  normal  to  the 
target  orbit  plane  (n)  can  then  be  defined  as 


n  = 


|r|  I^I 


(160) 


The  current  position  (^p)  and  velocity  (Rp)  of  the  payload  satellite 
also  can  be  determined  through  propagation  of  the  orbital  elements  specified 

at  epoch  for  that  satellite.  The  projection  (Y  )  of  the  satellite  position 

A  ^  ^ 

vector  (Rp)  into  the  plane  of  the  target  orbit  is  defined  as 

Yp  =  Rp  -  (tp  •  n)n  (161) 
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The  angular  separation  of  R  and  Y^,  measured  from  R,  can  then  be 


calculated 


X  = 


cos 


-1 


-»  -* 

Yp  •  R 

\%\  |R| 


^  ■¥ 

R  X 


|5pl 


•  n 


(162) 


Using  X,  the  true  anomaly  of  the  projection  Y  with  respect  to  the 
perigee  of  the  reference  orbit  (fp)  is  then  calculated 


fp  =  f  +  X 


(163) 


where  f  is  the  true  anomaly  of  target  (calculated  from  R  and  R).  The  radial 
distance  to  the  intersection  of  Yp  and  the  reference  orbit  trajectory  is 
given  by 


j.  _  a(l  -  e  )  ^ 

P  ~  1  +  e  cos  f 

P 

where 

a  =  semi-major  axis  of  the  reference  orbit 
e  =  eccentricity  of  the  reference  orbit 


(16^) 


To  evaluate  the  function  qp(t),  it  is  necessary  to  use  incomplete 
elliptic  integrals  of  the  second  type,  which  are  defined  in  terms  of  a  central 
angle  rather  than  a  focal  angle  (see  Fig.  49).  The  central  angle  |3 
corresponding  to  the  true  anomaly  of  the  reference  f  is  calculated  in  three 
steps : 

=  (ae)^  +  1r|^  +  2aejR|  cos  f 

ae  .  , 

sin  Y  =  sin  f 

3  =  f  -  Y  (165) 


( 
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Figure  49.  Relation  Between  Central  Angle  and  Focal  Angle 

A  similar  procedure  is  used  to  find  the  central  angle  associated  with  the 
true  anomaly  of  Yp 


2  2  2 

Q  =  (ae)  +  (rp)  +  2ae  Tp  cos  fp 
sin  Y  =  “q  sin  f., 

(ip  =  fp  -  Y  (166) 

The  arc  displacement  AL  between  the  center  of  the  debris  cloud  and  the 
projection  of  the  satellite  position  vector  onto  the  reference  orbit  plane  can 
then  be  expressed  as 


AL  =  a[E((5,  e)  -  E(|5p,  e)]  (167) 

where  the  function  E(3,  e)  is  defined  as  the  incomplete  elliptic  integral  of 
the  second  type 


E((l,  e) 


2  .  2 
e  sin 


d  ^ 
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For  small  values  of  eccentricity.  Reference  22  provides  a  good 
mate  method  for  evaluating  the  integral 


E(P,  e)  =  |(1  +  M)  ^^^[^215  +  I  sin  23 

„2 

-  (23  +  sin  23  cos  23) 

+  i,8  (2  sin  23  +  sin  23  cos^  23) 


+  0(e8) 


where 


M  = 


2 

e 


2  - 


2 

e 


The  function  qjCt)  (see  Fig.  50)  may  then  be  evaluated  as 


qj^(t)  =  |Al|  -  Lt 


Figure  50.  Determination  of  Satellite  Position 
Relative  to  Planar  Debris  Cloud 


approxi- 


(168) 


(169) 
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If  the  angle  y  is  positive,  then  the  payload  satellite  is  closest  to  the 

leading-edge  boundary  of  the  cloud  and  L  should  be  set  equal  to  the  leading- 
•  • 
edge  growth  rate,  .  If  the  angle  y  is  negative,  then  L  should  be  set  equal 

to  L^,  the  trailing-edge  growth  rate. 

Whenever  the  function  q^^  is  less  than  zero,  the  projection  of  the  payload 
satellite  position  vector  is  within  the  along-track  limits  of  the  debris 
cloud.  Once  the  torus  has  closed  (t  >  q^  is  always  negative. 

The  second  necessary  condition  for  the  payload  satellite  to  be  inside  the 
debris  cloud  is  that  the  function  q^  be  nonpositive.  The  previously  calcu¬ 
lated  quantities  Rp,  ?p,  ?p,  y,  and  rp  will  be  used  in  the  evaluation  of  q^. 

The  angular  distance  (9p)  from  Yp  to  the  point  where  the  target 
spacecraft  broke  up  is  calculated  by 

Gp  =  fp  -  f*  (170) 

where  f*  is  the  true  anomaly  of  the  point  where  the  spacecraft  broke  up. 

A  vector  (Rq)  from  the  intersection  of  ?p  with  the  target  reference  orbit 
trajectory  to  the  position  of  the  payload  satellite  may  be  defined  as 


-> 


The  elevation  angle  (^)  between  and  Yp  may  also  be  calculated 


cos  ^ 


(172) 
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Using  Eqs.  (145)  and  (146),  one  may  calculate  the  axes  of  the  elliptical  cross 
section  a(t,  0p)  and  b(t,  0p).  The  distance  (rg)  from  the  reference  orbit 
trajectory  along  to  the  boundary  of  the  cross  section  becomes 

2  u2r,  2  ^2,  -2  2  ,  ,, 

rg  =  b  [1  -  (a  -  b  )a  cos  (173) 

The  function  q2(t)  in  Figure  51  may  then  be  evaluated  as 


-  l^j)l  Tg 


(174) 


Figure  51.  Determination  of  Satellite  Position 

Relative  to  Debris  Cloud  Cross  Section 

8.5  DETERMINATION  OF  COLLISION  PROBABILITY 


Once  the  set  of  time  intervals  has  been  determined  during  which  the 
satellite  is  within  the  debris  cloud,  the  probability  of  collision  associated 
with  each  transit  can  be  calculated.  This  probability  is  a  function  of  satel¬ 
lite  cross-sectional  area,  time  spent  in  the  cloud,  relative  velocity  of  the 
satellite  with  respect  to  the  debris  particles,  average  cloud  density,  and  the 
total  number  of  debris  particles. 
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The  cumulative  probability  of  collision  (P^)  for  a  satellite  which 
travels  through  the  i-th  debris  cloud  a  number  of  times  can  be  defined  as 


M 

P.  =  n  (1  -  P.)  (175) 

"  j=l  ^ 


where  P^  is  the  probability  of  collision  associated  with  the  j-th  passage 
through  the  cloud,  and  M  is  the  total  number  of  passes  through  the  cloud. 


For  improved  computational  accuracy,  Eq.  (175)  may  be  evaluated  as 

P.  =  1  -  EXP 
1 

If  an  aggregate  of  debris  clouds  is  being  used  to  represent  a  larger 
cloud,  then  the  procedures  described  previously  must  be  applied  separately  to 
each  sub-cloud.  The  resulting  sets  of  time  intervals  are  then  used  to  calcu¬ 
late  a  cumulative  probability  of  collision  associated  with  each  sub-cloud 
[Eq.  (176)].  The  overall  probability  of  collision  corresponding  to  the  entire 
debris  cloud  (P)  may  then  be  calculated  as 


M 

I 


ln(l  -  P.) 


(176) 


P 


1  -  EXP 


N 

I  ln(l  -  P.) 
i=l  ^ 


(177) 


where  P^  is  the  cumulative  probability  of  collision  associated  with  the  i-th 
sub-cloud,  and  N  is  the  total  number  of  sub-clouds. 

8.6  EQUATIONS  FOR  DETERMINING  COLLISION  PROBABILITY  FOR 
PASSAGE  THROUGH  THE  DEBRIS  CLOUD 


This  section  describes  the  equations  used  to  calculate  the  probability 
of  collision  associated  with  a  particular  passage  of  a  satellite  through  a 
debris  cloud. 
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Assume  that  a  satellite  enters  a  debris  cloud  at  time  t^  and  exits  the 
cloud  at  time  tj^.  The  position  vector  of  the  satellite  at  these  two  times 
[^p(t^),  ^p(  tj^)]  can  be  determined  using  methods  previously  described.  The 
difference  in  true  anomaly  swept  out  by  passage  through  the  cloud  may  then  be 
calculated  as 


a  =  cos 


-1 


^p(to)  • 


^p(to)  1  I  l^p(t^) 


+  IK  INT[(t„-  t  )(T)~^] 
N  O 

where  T  is  the  period  of  the  reference  orbit. 


(178) 


The  number  of  sub-intervals  into  which  a  will  be  divided  may  then  be 
calculated  as 


N  =  1  +  INT 


/lO  \  -1 

®U8o)’^ 


(179) 


An  interval  of  10*  was  chosen  empirically. 


The  actual  distance  travelled  by  the  satellite  while  inside  the  cloud 
(dg)  may  be  approximated  as 


dg  =  i  a  |5p(t.)|  (180) 

where  t.  =  t  +  (i/N)  (t.,  -  t  ). 

10  No 

During  that  same  interval,  the  distance  travelled  by  the  cloud  (d-) 
is  approximately  equal  to 


-  '‘n- 


)  s  I 


1=0 


(181) 
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The  quantity  is  the  average  velocity  of  debris  particles  near 

the  satellite  at  time  and  is  given  by 


Vp(t^)  =  \J\iillr  -  1/a) 


(182) 


where  a  is  the  semi-major  axis  of  the  reference  orbit,  p  is  the  gravitational 

constant,  and  r  is  calculated  from  Eq.  (164). 

P 

The  relative  distance  traveled  by  the  satellites  with  respect  to  the 

debris  cloud  (d  , )  becomes 
rel 


d  ,  =  d„  -  cos(Ai)d„  (183) 

rei  b  L 

where  Ai  is  the  difference  in  inclination  between  the  satellite  and  target 
orbit  planes. 

The  effective  volume  swept  out  by  the  satellite  moving  through  the 
cloud  (VOLp)  may  then  be  calculated  as 

VOLp  =  Ap|d^^J  (184) 

where  Ap  is  the  cross-sectional  area  of  the  satellite.  The  average  volume 
of  the  debris  cloud  associated  with  the  region  of  passage  (VOL^^)  can  be 
found  using  either  method  previously  described.  These  volumes  are  used  to 
compute  the  probability  of  collision  associated  with  this  partial  passage 
through  the  cloud: 


P 


1  -  EXP 


Np  In  I 


where  Np  is  the  total  number  of  debris  particles  in  the  cloud. 


(185) 
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9.  EFFECTS  OF  ECCENTRICITY  ON  THE  VOLUME  OF  A  DEBRIS  CLOUD 

9.1  INTRODUCTION 

During  evaluation  of  the  collision  hazard  posed  by  the  debris  from  an 
orbital  breakup  of  a  spacecraft,  it  is  necessary  to  define  a  region  that  the 
debris  would  occupy.  This  region  is  a  time-varying  volume,  with  each  particle 
occupying  a  different  orbit  and  all  spreading  out  relative  to  each  other. 

The  first  approximation  to  this  debris  cloud  volume  assumes  that,  prior 
to  breaking  up,  the  spacecraft  in  question  is  in  a  circular  orbit  about  the 
earth.  Additionally,  the  linearized  equations  of  motion  are  used,  remaining 
valid  for  small  changes  in  relative  velocities  and  over  a  small  number  of 
orbits . 


In  this  section,  small  values  of  eccentricity  are  added  into  the 
equations  of  motion.  By  using  a  differential  correction  process,  similar  to 
the  work  in  Reference  23,  one  finds  new  functions  involving  the  circular 
solution  and  the  changes  in  that  solution  due  to  eccentricity.  Additionally, 
the  solution  process  allows  for  the  orbital  breakup  to  occur  at  any  time 
throughout  the  elliptical  orbit. 

This  continuing  process  to  refine  the  volume  model  leads  to  a  greater 
understanding  of  what  is  really  happening,  and  allows  for  a  more  accurate 
assessment  of  the  collision  hazard  posed  by  the  debris. 

9.2  ANALYSIS 


The  solution  process  involves  a  number  of  fundamental  techniques,  as 
well  as  some  new  ideas  first  introduced  for  this  problem.  Initially,  the 
equations  of  motion  from  Newton's  laws  are  formulated;  from  there,  the 
solution  is  perturbed  by  adding  eccentricity.  The  equations  are  linearized, 
and  the  state  transition  matrix  relating  position  and  initial  relative 
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velocity  is  found.  This  procedure  was  used  by  Anthony  and  Sasaki  (Ref.  24). 
Finally,  using  a  method  developed  specifically  for  this  problem,  one  finds  the 
volume  of  the  debris  cloud. 

Given  an  inertial  frame  (I)  with  origin  at  the  center  of  the  earth,  and 
a  rotating  frame  (R),  with  Y  being  in  the  radial  direction  and  Z  being  in  the 
out-of-plane  direction  as  seen  in  Figure  52,  the  position  vector  of  the  target 
and  i-th  piece  of  debris,  respectively,  is 

?  *  r  J  (186 


R^  =  XI  +  (Y  +  r)  J  +  ZK 


while  the  vector  between  the  target  and  the  debris  is 


R  =  R.  -  r  =  XI  +  YJ  +  ZK 
1 


VEHICLE 
WITHOUT 
BREAKUP -r 


i-th  DEBRIS  PARTICLE 
(after  breakup) 


INERTIAL 

REFERENCE 

FRAME 


Figure  52.  Coordinate  Frame 


The  differential  equations  in  the  inertial  frame  for  this  inverse  square  field 
are 


dt 


(189) 


The  left  side  of  Eq.  (189)  is  obtained  by  taking  the  derivative  of  Eq.  (188) 
in  the  inertial  frame 


dt  ~  dt 


■«-  (I)  X  R 


(190) 


where 


and 


_(R)  « 

dR^  ^  • 

-  —  -  =  XI  +  YJ  +  ZK 
dt 


«  =  §K 


(191) 


(192) 


The  second  derivative  of  Eq.  (190)  in  the  inertial  frzune  is 


d^R^^^  d 


dt 


dt 


dt 


+  —  (w  X  R) 


(193) 


When  the  right-hand  side  of  Eq.  (189)  is  expanded  and  the  individual  compo¬ 
nents,  are  equated,  the  equations  of  motion  become 


X  -  Y9  -  2y§  -  + 


M. 


IX^  +  (Y  +  r)^  + 


Y  +  xe  +  2X5  -  §^Y  +  — z - ^ 


[X^  +  (Y  +  r)^  + 


=  0 


Z  + 


uZ 


[X^  +  (Y  +  r)^  + 


=  0 


(194a) 


(194b) 


(194c) 
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Nondlmensional  terms  are  introduced  such  that 


X 


y 


z 


p 


X 

a 

Y 

a 

Z 

a 

r 

a 


where  a  is  the  semi-major  axis  of  the  reference  orbit,  so  Eqs. 


x"  -  ye"  -  2y'0'  -  0'^x  + 


,  ,  v2  2,3/2 

[x  +  (y  +  p)  +  z  J 


y"  +  xe"  +  2x'e’  -  0’^y  + 


i3L-l-0l 


,  2  t  ,2  2,3/2 

[x  +  (y  +  p)  +  z  ] 


X  2  ,  ,22,  3/2 

[x  +  (y  +  p)  +  z  ] 


=  0 


(195) 

(196) 

(197) 

(198) 

(199) 

(194)  reduce  to 
=  0  (200a) 

=  0  (200b) 

(200c) 


where  '  (prime)  represents  differentiation  with  respect  to  t.  Assuming  that 
the  distance  between  a  debris  particle  and  the  target  is  much  less  than  the 
semi-major  axis,  one  finds  that  Eqs.  (200)  further  reduce  to 


x"  -  y0"  -  2y'0’  + 


4 

P 


y"  +  x0"  +  2x'0’ 


(x^  -  2y^  +  z^)  =  0 
2P 


z 


It 


0 


(201a) 


(201b) 


(201c) 
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For  a  slightly  elliptical  orbit*  the  series  for  0  and  R  are 
/_y_\l/2  r,  ^  w  5  2 


0  = 


[1  +  2e  cos  M  +  ^  e^cos  2M  +  ...] 


r  =  (— )  [1  -  e  cos  M  +  §  (1  -  cos  2M)  +  ...] 

CO  2, 


(202) 


where 


1/2 


(t  -  t^) 
P 


(  3) 

'  a  ' 

Nondimens ionally,  and  retaining  only  linear  terms,  Eq.  (202)  becomes 


(M  -  mean  anomaly) 

(tp  -  time  of  periapsis  passage) 


0'  =  1  +  2e  cos(t  -  Tp) 

p  =  1  -  e  cos(t  -  tp)  (203) 

For  the  circular  target  orbit,  the  Clohessy-Wiltshire  equations  take  the  form 


Xc" 

-  2yc  =  0 

(204a) 

yc" 

+  2xc  -  3yc  =  0 

(204b) 

2c" 

+  Zc  «  0 

(204c) 

Now,  assuming  a  new  solution  for  x,  y,  and  z,  in  the  form 


X  =  Xc  +  6x  (205a) 
y  =  yc  +  5y  (205b) 
z  =  Zc  +  5z  (205c) 


When  one  takes  appropriate  derivatives,  while  retaining  only  linear  terms,  the 
new  differential  equations,  in  terms  of  the  perturbations  on  x,  y,  and  z,  are 

8x"  -  25y'  =  e[(4y^  +  x^)cos(t  -  x^)  -  2y^sin(T  -  x^)]  (206a) 

6y"  +  28x'  -  38y  =  e[(10y  -  4x')cos(x  -  x  )  +  2x  sin(x  -  x  )]  (206b) 

c  c  pc  p 

8z"  +  8z  =  -  3ez  co8(x  -  x  )  (206c) 

c  p 
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Initially 


Sx(0)  =  5y(0)  =  Sz(0)  =  6x’(0)  =  6y’(0)  *  Sz'CO)  =  0 
Assuming  a  solution  of  Eqs,  (206)  in  the  form 

6x  =  etAp  +  Aj^sin  x  +  A2COS  x  +  A^sin  2t 

+  A^cos  2x  +  A^x  +  A^x  sin  x  +  A^x  cos  x]  (207a) 

Sy  =  e[BQ  +  B^sin  x  +  B^cos  x  +  B^sin  2x 

+  B^cos  2x  +  B^x  +  BgX  sin  x  B^x  cos  x]  (207b) 

5z  =  e[CQ  +  C^sin  x  +  C^cos  x  +  C^sin  2x 

+  C^cos2  X  +  CjX  +  CgX  sin  x  +  C^x  cos  x]  (207c) 

where 

1 

A«  =  -3x'  sin  X  -  o  y’  cos  x 
0  0  p  2  "^o  p 

A-  =  6y'  sin  t 

1  o  p 

A„  =  6x'  sin  T  +  2y*  cos  x 

2  o  p  o  p 

3 

A-  =  -  -ry'  sin  x  +  3x’  cos  x 

3  2"^  o  p  o  p 

3 

A,  =  -3x'  sin  X  -  r  y*  cos.  x 
h  o  p  2  "^0  p 

A5  =  -3yo  sin  Xp  -  3Xo  cos  Xp 

Ae  =  -3x^  sin  Xp 

A;  =  -3Xo  cos  Xp  (208) 

®0  =  -3yo  sin  Xp  -  4xo  cos  Xp 

®1  =  -Xq  sin  Xp  -  2yo  cos  Xp 

B2  =  ‘*yo  sin  Xp  +  2x^  cos  Xp 

B3  =  2xo  sin  Xp  +  y'o  cos  Xp  (209) 

(cont. ) 
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^4  =  -  yo  sin  tp  *  2Xo  cos  Xp 
B5  =  0 

Be  =  3xo  cos  Xp 
B7  =  -  3xq  sin  Xp 


C 

C 

C 

C 

C 


0 

1 

2 

3 

4 


'  sin  X 
o 

cos  X 


2z'  sin  X 
0  p 


2  o 


z  cos  X 


2  o 


z'  sin  X 


P 


€5=0 
Ce  =  0 
C7  =  0 


(209) 


(210) 


When  Eqs.  (208)  through  (210)  are  inserted  into  Eqs.  (207),  and  the  results 
placed  in  matrix  form 


5x 


=  [1  sin  X  cos  X  sin  2x 

cos  2x  X  sin  x  x  cos  x] 


[({>] 


-3e  sin  Xp 
0 

6e  sin  x 

P 

3e  cos  Xp 

-3e  sin  Xp 
-3e  cos  Xp 
-3e  sin  Xp 

-3e  cos  Xp 
[AE] 


1 

2®  nos 

■'P 

0 

r  1  ' 

Xq 

6e  sin 

■'P 

0 

' 

yo 

3 

2®  sin 

X 

P 

0 

1 

.  ^0 

3 

2®  sin 

^P 

0 

3 

■  2®  cos 

•'P 

0 

-3e  sin 

■'P 

0 

0 

0 

0 

0 
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(211a) 


In  a  more  compact  form,  Eqs.  (211)  become 


Similarly,  the  Clohessy-Wiltshire  equations  can  be  written  in  the  same  form 


where 


'  •> 

Xc 

= 

“[<l>] 

[ACf 

<■  ■< 

t 

Xo 

yc 

•  = 

[♦I 

[BC] 

yo 

2c 

V 

= 

[CC], 

1 

2o 

[AC] 


0  2  0 

4  0  0 

0-2  0 
0  0  0 

0  0  0 

3  0  0 

0  0  0 

0  0  0 


(213) 


(214a) 


[BC] 


-2  0  0 

0  10 

2  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 


(214b) 


[CC] 


0  0  0 
0  0  1 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0. 


(214c) 
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Combining  these  two  solutions 


f  ' 

6x 

W  [A] 

Xo 

•  _  < 

yc 

■  +  ' 

Sy 

'  = 

[♦]  [B] 

y; 

) 

.  ^c. 

. 

W  [C]_ 

^o 

1 

= 

[M*] 

1 

where 


[A]  =  [AC]  +  [AE] 


-3e  sin  Xp 
k 

6e  sin  Xp 
3e  cos  Xp 

-3e  sin  Xp 

-3  -  3e  cos  Xp 
-3e  sin  Xp 
_23e  cos  Xp 


2  -  2  e  cos  Xp 

6e  sin  Xp 

-2  +  2e  cos  Xp 
3 

-  2  e  sin  Xp 
3 

-  2  e  cos  Xp 

-  3e  sin  Xp 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


(215) 


(216a) 


[B]  =  [BC]  +  [BE] 


-2  -  4e  cos  Xp 
-e  sin  Xp 
2  +  2e  cos  Xp 
2e  sin  Xp 
2e  cos  Xp 
0 

-3e  cos  Xp 
j^3e  sin  Xp 


-3e  sin  Xp 
1  -  2e  cos  Xp 
4e  sin  Xp 
e  cos  Xp 
-e  sin  Xp 
0 
0 
0 


0 

0 

0 

0 

0 

0 

0 

0 


(216b) 
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0  0 

0  0 

0  0 

[C]  =  [CC]  +  [CE]  =0  0 

0  0 

0  0 

0  0 

0  0 


The  volume  of  a  debris  cloud 


-  -e  sin  Tp 
le  cos  Tn 


2e  sin  t. 


-e  cos  Tp 


-  -e  sin  Tp 


0 

0 

0 


V  =  ^  det 


1/2 


(216c) 


(217) 


This  differential  correction  process  is  very  useful  and  can  be  utilized  for 
many  other  types  of  perturbations,  simply  by  getting  the  perturbation  in  the 
form  of  Eq.  (211)  and  adding  the  8x3  matrices  together  with  the  linear 
solution  of  Eq.  (216). 

9.3  RESULTS 


The  effects  of  location  in  an  elliptic  orbit  are  shown  for  various 
values  of  eccentricity  from  0  to  0.25  in  Figures  53  through  60.  For  perigee 
and  apogee  breakups,  the  volume  profiles  are  very  similar  to  a  circular  case. 
However,  at  points  halfway  between  apogee  and  perigee  (x^  =  90*,  270®), 
there  is  a  distinct  deviation  for  different  eccentricities.  As  e  approaches 
0.25,  the  analysis  deteriorates  and  produces  doubtful  results. 
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CONCLUSIONS 


9.U 


As  shown  in  the  results  of  this  analysis,  the  addition  of  eccentricity 
to  the  initial  orbit  of  a  spacecraft  prior  to  orbital  breakup  has  a  distinct 
effect  on  the  evolution  of  its  debris  cloud.  For  very  small  changes  in  eccen¬ 
tricity  (from  e  =  0.01  to  e  ®  0.05),  subtle,  yet  distinctive,  characteris¬ 
tics  are  present.  At  apogee  and  perigee  breakups,  the  half-orbit  zero  points 
are  the  same  as  the  circular  case;  but  when  the  breakup  occurs  at  some  other 
time,  that  half-orbit  zero  point  moves  away  from  exactly  half  an  orbit. 

Although  not  obvious  for  the  e  =  0.01  case,  this  half-orbit  point  shift  becomes 

more  apparent  as  eccentricity  increases.  Additionally,  as  eccentricity 

2 

reaches  0.25,  the  assumption  that  e  terms  are  dropped  becomes  invalid  and 
produces  doubtful  results. 


i 


I 


I 


Volum*  (non-dimensional)  >13  Volume  (non-dimensional) 


Volume  vs.  Number  of  Orbits 

T,  =  180 


Lgure  57. 


Volume  versus  Time;  x 
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180*  (apogee) 


Volume  vs.  Number  of  Orbits 

T,  =  225 


Figure  58. 


Volume  versus  Time;  x 

P 


225  • 
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Volume  vs.  Number  of  Orbits 


=  270 


Figure  59.  Volume  versus  Time;  =  270® 


Volume  VS.  Number  of  Orbits 

T,  =  315 


Figure  60.  Volume  versus  Time;  x  *  315® 


APPENDIX  A 

DEBRIS  CLOUD  VOLUME  AS  A  FUNCTION  OF  TIME 

The  determinant  of  M  in  Eq.  (4)  is  the  volume  of  a  parallelepiped  P  in 
three-dimensional  space,  when  the  edges  of  P  are  obtained  from  the  rows  of  M  as 
illustrated  in  Figure  A-1  (Ref.  A-1).  To  show  this,  suppose  that  all  angles  in 
Figure  A-1  are  right  angles  or  that  the  edges  of  P  are  perpendicular.  Then  the 
volume  is 


VOLUME  =  1^  1^  I3  (A-1) 

Since  the  edges  of  P  are  the  rows  of  M  (or  columns  with  a  different  P 
but  the  same  volume)  which  are  mutually  perpendicular,  the  orthogonality  of 
the  matrix  M  yields 


T 

M  M  = 


(A-2) 


2  2  2 

The  Ij^  ,  1^  »  I3  are  the  squares  of  the  lengths  of  the  rows,  and  the 
zeroes  off  the  diagonal  are  the  result  of  orthogbnality  of  the  rows.  Taking 
determinants 


det  [M^M  ]  *  [det  M]  [det  M^J 
*  [det 


S 


1 


2 

3 


=  (VOLUME)^ 


(A-3) 


If  the  region  is  not  rectangular,  M  is  not  an  orthogonal  matrix  and  the 
volume  is  no  longer  the  product  of  edge  length.  However,  a  transformation  can 
be  found  which  diagonalizes  M;  this  is  analogous  to  finding  the  "principal"  or 
orthogonal  axes  for  the  system.  The  corresponding  elements  of  the  diagonal 


A-1 


1 


y 


Figure  A-1.  Volume  of  Determinant  M 

matrix  are  the  required  lengths  of  P.  They  can  be  determined  as  follows: 
find  the  roots  of  M  by  solving  the  characteristic  equation 

|M  -  Xlj=  0  (A-^) 


where  I  is  the  unit  matrix. 

For  the  case  in  Eq.  (5),  Section  2,  the  characteristic  equation  is  of  the 

form 


[(a^^  -  X)  (a22  -  X)  -  a^^]  ia^j  -  X)  =  0  (A-5) 

The  roots  of  the  characteristic  equation  are 


\  “  ®33 

(a^l  +  a22)  ±  +  322)^  -  Ma^^  a22 

X2  3  -  - 2 


^2  ^21^ 


1/2 


(A-6) 


A-2 


where  the  X„  -  roots  are  obtained  from  the  bracketed  term  in  Eq.  (A-5).  The 
product  of  the  three  roots  X^,  X^,  X^  represents  the  volume  of  the  cloud.  Since 
the  volume  is  a  positive  scalar  quantity,  the  product  of  the  roots  must  be  a 
real  positive  value.  The  necessary  and  sufficient  conditions  are  that 


X 


1 


(A-7) 


and 


X^  •  X^  >  0  (A-8) 

The  second  condition  can  be  satisfied  if  X^  and  X^  are  complex  .conjugates 
(since  the  product  of  two  conjugate  complex  systems  is  positive)  or  both 
^2’  ^3  positive  real  numbers.  Since  the  roots  represent  the  sides 

of  P,  only  the  condition 


(A-9) 


needs  to  be  satisfied.  This  requires  that  the  terms  in  the  radical  satisfy 
the  condition 


(a 


11 


>  4(a 


11  22 


-  a^2  ^21) 


(A-10) 


or  that 


®11  ®22 


®12  ®21 


>  0 


(A-11) 


In  terms  of  the  matrix  [Eq.  (5),  Section  2],  it  is  necessary  and 
sufficient  that 


ad  +  b  >  0 


(A-12) 


A-3 


Physically,  this  means  that  the  area  of  the  cloud  in  the  x,y  plane  is 
always  finite  and  positive  except  at  9  =  0,  360®,  et  cetera,  when  it  is 
zero.  The  normalized  volume  of  the  cloud  is  of  the  form 


VOLUME  =  det  [M] 


(A-13) 


where 


det[M]  =  1 (a, .  a,-  -  a 


11  22  12  “21 


=  I  (ad  +  b  )  •  (dj 


=  ft  Q,  Q, 
12  3 


)|  •  1^331 


Equation  (A-13)  is  the  volume  of  a  spheroid  with  the  serai-major  axes 
ft^,  ft^,  ^3  resulting  from  the  application  of  a  unit  velocity  impulse  to 
particles  with  uniform  (spherical)  distribution  in  the  x,y,z  reference  frame. 
By  analogy,  the  area  of  the  cloud  in  the  x,y  plane  can  be  expressed  as 


AREA  =  irft.ft. 
x,y  1  2 


=  ir  ad  +  b 


(A-14) 


This  is  for  a  circular  distribution  of  the  \init  velocity  impulse  in  the  x,y 
plane. 


APPENDrX  B 


UNIFORMLY  DISTRIBUTING  POINTS  ONTO  A  SPHERE 


B.l  INTRODUCTION 

In  the  development  of  a  collision  breakup  model,  a  basic  problem  is  to 
determine  the  resulting  fragments'  directions.  In  our  model,  the  fragment 
directions  are  chosen  randomly  from  a  uniform  distribution  of  directions. 

This  appendix  describes  the  method  for  determining  the  uniform  distribution  of 
directions.  These  directions  are  equivalent  to  points  projected  onto  a 
sphere.  Uniform  distribution  then  means  that  all  adjacent  point-  on  the 
sphere  are  equidistant  from  one  another.  The  method  can  be  applied  to  both 
collision  and  explosion  models. 

Following  a  description  of  the  method  to  uniformly  spread  points  onto  a 
sphere,  two  examples  of  its  use  will  be  presented  using  a  satellite  explosion. 
A  satellite  in  an  orbit  which  has  a  semi-major  axis  of  A444  nmi  is  assumed  to 
explode  uniformly,  imparting  a  velocity  of  1000  ft/sec  on  all  the  fragments. 

In  the  first  example,  the  satellite  is  in  a  circular  orbit;  in  the  second,  it 
is  in  an  elliptical  orbit.  The  resulting  orbits  for  each  fragment  are 
computed  using  the  On-Line  Orbital  Mechanics  (OLOM)  program. 

B.2  ANALYSIS 

The  method  for  distributing  points  uniformly  onto  a  sphere  is  the  same 
as  that  used  for  constructing  geodesic  domes.  The  method  starts  with  the 
fifth  regular  polyhedron  (platonic  solid):  an  icosahedron.  Convex  polyhedra 
are  considered  regular  if  they  have  faces  that  are  regular  polygons  and  all 
the  polyhedral  angles  are  congruent. 

An  icosahedron  has  12  vertices  and  20  faces,  which  are  equilateral 
triangles  (Fig.  B-la).  The  distance  from  any  vertex  to  the  geometric  center 
is  the  same.  Hence,  if  a  sphere  were  circumscribed  about  an  icosahedron,  its 
geometric  centers  would  coincide  and  the  12  vertices  would  all  lie  on  the 
surface  of  the  sphere.  Since  the  distances  between  two  adjoining  vertices  of 
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a  regular  polyhedron  are  equal,  it  follows  that  the  arc  lengths  on  the  surface 
of  the  sphere,  joining  any  two  adjacent  vertices,  also  are  equal.  Thus, 
simply  by  circumscribing  a  sphere  about  an  icosahedron,  we  have  found  12 
points  uniformly  distributed  on  a  sphere. 

In  order  to  distribute  more  points  about  the  sphere,  first  establish  the 
origin  of  a  three-dimensional,  rectangular  coordinate  system  at  the  geometric 
center  of  the  polyhedron  as  shown  in  Figure  B-lb  (which  is  Fig.  6  in  Ref.  B-1). 
This  orientation  is  for  geometrical  calculations.  Table  B-1  lists  the  coordi¬ 
nates  of  the  12  vertices  of  the  principal  polyhedral  triangles  (PPT).  A  PPT 
is  any  one  of  the  equal  equilateral  triangles  which  forms  the  face  of  the 
regular  polyhedron  (Ref.  B-1). 

Table  B-1.  Coordinates  of  the  PPT's  Vertices  of  an  Icosahedron 


Vertices; 


(0,  +  7x75^^^,  +  1/5^^^  /x) 


where 


(+  A,  0,  + 


(+  X,  0) 


X  =  =  1.61803 


Pj  »  (0,  /t/5^^^,  1/5^^^  Vx)  =  (0,  0.85065081,  0.52573111) 


P,  =  (1/5^^^  \/x,  0,  %/x/5^^^)  =  (0.52573111,  0,  0. 


85065081) 


P»  =  (V^/5^^^,  1/5^''^V^,  0)  =  (0.85065081,  0.52573111,  0) 
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Figure  B-la.  Icosahedron 


Figure  B-lb.  Three-Dimensional  Coordinate  System  at 
the  Geometric  Center  of  the  Polyhedron 


The  next  step  is  to  subdivide  the  PPT  into  frequency  N  as  shown  in 
Figure  B-2a.  The  frequency  is  the  number  of  parts  or  segments  into  which  a 
principal  side  is  subdivided.  A  principal  side  is  any  one  of  the  three  sides 
of  the  PPT  (Ref.  B-1).  Figure  B-2b  shows  how  each  subdivision  point  is 
connected  by  line  segments  parallel  to  their  respective  sides.  A  grid  of 
equilateral  subtriangles  is  generated. 


The  coordinates  of  the  equilateral  subtriangle  vertices  (Xjj, 


Z^.)  are  obtained  from  Eq.  (B-1) 

X  J 


X  -  X  X  -  X 

=  ^1 " '  "  -  N— 


+ 


Z 


IJ 


+ 


(B-1) 


where  N  is  the  frequency,  (Xj^,  Yj^,  Z^)  (.^2*  ^2’  ^2^  (^3*  ^3*  ^3^ 
coordinates  of  the  PPT  vertices,  and  I  and  J  are  integers  such  that  0  <  J 
<  I  <  N  (Ref.  B-2).  The  I  and  J  have  unique  values  for  each  vertex;  they 
label  each  vertex  as  shown  in  Figure  B-3  (which  is  Fig.  10  in  Ref.  B-3). 


For  the  purposes  of  this  study,  it  is  necessary  only  to  calculate  each 
vertex  direction  from  the  origin  and  not  the  vertex  distance  from  the  origin. 
This  is  calculated  by  doing  a  coordinate  transformation  from  rectangular  to 
spherical  coordinates  for  each  subtriangle  vertex  of  each  PPT.  The  resulting 
spherical  coordinates  of  each  vertex  direction  (0,  ^)  represent  the  direc¬ 
tions  in  which  Av's  are  applied  for  each  particle  of  the  exploded  satellite. 
These  coordinates,  along  with  specified  orbital  parameters  and  a  specified 
Av,  are  input  to  OLOM.  Orbital  parameters  are  computed  for  each  particle 
from  the  exploded  satellite. 
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Figure  B-2a.  Subdivided  PPT  into  Frequency  N 


Figure  B-2b.  Grid  of  Equilateral  Subtriangles  in  PPT 


B.3  EXAMPLE  1 


In  this  example,  the  satellite  is  assigned  to  be  in  a  circular  orbit  with 
a  1000-nmi  (a  =  4A44  nmi)  altitude.  The  inclination  is  28.5®,  and  the 
satellite  is  at  its  southernmost  point  when  it  explodes.  The  Av  applied  to 
each  particle  from  the  explosion  is  1000  ft/sec.  A  frequency  of  N  equal  to  7 
is  used  to  generate  the  particles. 

One  interesting  result  is  obtained  by  graphing  each  particle's  resulting 
apogee  and  perigee  against  its  period  (Fig.  B-A).  This  result  was  origin¬ 
ally  demonstrated  by  John  Gabbard  (Ref.  B-3).  Figure  B-4  has  two  wings  that 
meet  at  a  point.  The  upper  wing  shows  the  fragment  apogees  plotted  against 
their  periods,  and  the  lower  wing  shows  the  perigees  versus  the  periods.  Two 
points  are  plotted  for  each  fragment.  Note,  also,  the  uniform  distribution  of 
particles  in  this  figure. 

To  calculate  the  number  of  particles  generated  by  this  method,  recall 
that  an  icosahedron  has  12  vertices,  20  faces,  and  30  edges.  If  a  PPT  is 
subdivided  into  frequency  N,  then  the  number  of  particles  (NP)  generated  is 


N-2 

NP  =  12  +  20  X  i  +  30(N  -  1)  (B-2) 

i=l 


5 

Thus,  for  N  =  7,  NP  =  12  +  20  ^  i  +  (30  x  6)  =  492  particles  are  generated. 

i=l 

However,  since  Figure  B-4  plots  apogee  and  perigee  versus  period  for  each 
particle  from  the  explosion.  Figure  B-4  has  984  points  plotted. 

In  order  to  more  fully  understand  the  relation  between  the  location  of  a 
point  on  the  curve  and  the  actual  direction  of  the  Av  that  generated  that 
point,  Av's  of  1000  ft/sec  were  applied  in  10*  increments  in  three  orthogonal 
planes.  The  planes  are  the  yaw,  pitch,  and  roll  planes  shown  in  Figure  B-5. 
The  yaw  plane  is  defined  by  the  velocity  vector  ^  and  the  momentum  vector 
r  X  V.  The  pitch  plane  is  defined  by  the  radius  r  and  velocity  V  vectors. 

The  roll  plane  is  perpendicular  to  the  other  two  planes. 
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/ vvv , ,  /vv.\ 

(N.  0)  (N,  1)  (N,  2)  (N.  3)  (N,  N-3)  (N,  N-2)  (N.  N-1)  (N,  N) 


Figure  B-3.  Breakdown  Numbering 


PERIOD  (min) 


Figure  B-4.  Fragment  Apogees  and  Perigees  versus  Their 

Periods  from  a  Satellite  Exploding  Uniformly 
(A  =  4A4A  nmi;  ECC  *  0;  NF  =  492) 
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Figure  B-5  also  shows  that  the  0*  direction  in  both  the  yaw  and  pitch 
planes  is  in  the  V  direction.  The  0“  direction  for  the  roll  plane  is  in  the 
-r  X  V  direction.  For  both  the  roll  and  pitch  planes,  the  90“  direction  is  the 
r  direction,  and  for  the  yaw  plane  it  is  the  -r  x  V  direction. 

Figure  B-6  shows  the  resulting  apogees  and  perigees  from  a  1000  ft/sec 
Av  applied  in  10“  increments  in  the  three  orthogonal  planes  shown  in  Figure 
B-5.  Notice  first  that  the  results  in  both  the  yaw  and  pitch  planes  are 
symmetric  about  180“;  the  results  are  the  same,  for  example,  whether  the  angle 
is  130“  or  -130“.  The  roll  plane  is  S3rmmetric  about  90“  and  180“;  the  results 
are  the  same,  for  example,  whether  the  angle  is  50“,  -50“,  130“,  or  -130“. 


Also  note  that,  as  expected,  the  Av  applied  in  the  V  direction  (yaw  and 
pitch  direction  is  0“)  yields  the  highest  energy  orbit.  When  the  Av  is 
applied  solely  in  the  roll  plane,  the  resulting  orbit  period  is  nearly  that  of 
the  parent  body. 

The  final  point  is  made  by  comparing  Figures  B-4  and  B-6.  The  outline  in 
Figure  B-4  that  contains  all  points  is  the  same  outline  in  Figure  B-6  that  was 
produced  from  a  Av  applied  in  the  yaw  and  pitch  planes.  The  thickness  of 
the  outline  is  determined  in  the  roll  plane.  Thus,  an  idea  of  the  pitch,  yaw, 
and  roll  angles  of  the  applied  Av  can  be  determined  for  any  specific  point 
in  Figure  B-4. 

B.4  EXAMPLE  2 

The  assumptions  in  this  example  are  the  same  as  those  in  the  first,  except 
that  this  orbit  has  an  eccentricity  equal  to  0.05.  The  semi-major  axis  is 
4444  nmi,  the  inclination  is  28.5*.  The  Av  applied  to  each  particle  from 
the  explosion  is  1000  ft/sec.  Finally,  the  frequency  used  to  distribute  the 
particles  is  7. 
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Figure  B-5.  Angles  in  Three  Orthogonal  Planes 
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B-6.  Apogees  and  Perigees  versus  Period  for  Av  =  1000  ft/sec 
Applied  in  Three  Orthogonal  Planes  to  an  Object  in 
1000  nmi  Circular  Orbit 


Figures  B-7  and  B-8  show  the  resulting  apogees  and  perigees  of  the 
exploded  fragments  versus  their  periods  at  two  positions  in  the  orbit. 

Figure  B-7  shows  the  resulting  altitudes  if  the  explosion  occurs  at  perigee 
(h  =  778  nmi),  while  Figure  B-8  shows  the  altitudes  if  the  explosion  occurs  at 
apogee  (h  =  1222  nmi). 

An  interesting  observation  can  be  made  by  comparing  the  two  figures. 
Figure  B-7  shows  that  the  altitude  range  of  the  fragment  perigees  (resulting 
from  the  explosion  at  perigee)  varies  from  490  to  778  nmi.  Figure  B-8  shows 
that  the  altitude  range  of  fragment  perigees  (resulting  from  the  explosion  at 
apogee)  varies  from  100  to  1222  nmi.  The  median  perigee  altitude  for  both 
cases  is  near  650  nmi.  The  observation  to  be  made  is  that  the  perigee 
altitudes  go  much  lower  in  Figure  B-8  than  in  Figure  B-7.  Thus,  some  of  the 
particles  in  Figure  B-8  will  reenter  the  atmosphere  relatively  quickly, 
whereas  no  particles  will  reenter  quickly  in  Figure  B-7.  Therefore,  from  a 
debris  hazards  standpoint,  if  a  satellite  is  to  be  exploded  in  orbit,  a  case 
can  be  made  for  exploding  it  at  apogee  rather  than  perigee. 

B.5  SUMMARY 

This  appendix  presents  a  geometric  method  for  uniformly  distributing 
points  onto  a  sphere.  Two  applications  of  its  use  are  presented  using  a 
satellite  explosion.  The  first  application  is  a  satellite  explosion  in  a 
circular  orbit,  and  the  second  is  an  explosion  in  an  elliptical  orbit.  Graphs 
are  presented  showing  the  apogees  and  perigees  of  the  resulting  fragments 
versus  their  periods. 
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Figure  B-7.  Fragment  Apogees  and  Perigees  versus  Their  Periods 
from  a  Satellite  Exploding  Uniformly  at  its 
Perigee  (A  =  Uhhh  nmi;  ECC  =  0.05;  NF  =  492) 
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Figure  B-8.  Fragment  Apogees  and  Perigees  versus  Their  Periods 
from  a  Satellite  Exploding  Uniformly  at  its 
Apogee  (A  =  4444  nmi;  ECC  =0.05;  NF  =  492) 
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